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19.  Abstract  (continued) 

within  of  observations  by  a  seismic  network,  but  will  be  only  partially  successful 
in  reducing  the  effects  of  near-source  structure.-; 


1  Summary 
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1.1  OBJECTIVES 


Observations  of  teleseismic  P  waves  at  large  aperture  arrays  have  found  amplitude 


fluctuations  across  the  array  as  large  as  those  observed  over  worldwide  networks.  The 


amplitude  fluctuations  and  their  correlation  with  travel  time  variations  are  consistent 


with  the  focusing/defocusing  effects  of  three-dimensional  velocity  structure  beneath  the 


array,  rather  than  the  effects  of  intrinsic  attenuation.  A  reciprocal  effect  on  amplitudes 


is  to  be  expected  for  variations  in  the  location  of  underground  nuclear  tests  within  a  test 


site.  An  objective  of  this  project  is  to  predict  these  amplitude  variations  using  known  3- 


D  structure  within  a  test  site.  Fast,  asymptotically  approximate  methods  of  modeling 


the  body  wavefield  are  used  to  predict  these  amplitude  variations  by  forward  modeling 


in  structures  determined  from  block  3-D  inversions  of  travel  times,  as  well  as  more 


detailed  models  determined  from  local  geophysical  surveys  at  the  Nevada  test  site.  The 


results  of  these  studies  can  be  used  to  determine  the  fundamental  scale  lengths 


required  to  produce  a  correctable  amplitude  anomaly  of  a  given  size. 


1.2  RESULTS 


1.2.1  NTS 


A  known  three-dimensional  structure  beneath  Pahute  Mesa,  Nevada  Test  Site 


(Taylor,  1983)  can  account  for  many  of  the  features  in  the  azimuthal  amplitude  pattern - 


of  teleseismic  P  waves  from  Pahute  underground  tests.  This  model  can  be  used  to 


correct  for  focusing  and  defocusing  effects  of  the  structure  beneath  Pahute  Mesa 
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accounting  for  factors  of  three  in  amplitude  fluctuation  and  for  0.6  sec.  in  travel  time 
fluctuation.  The  reduction  in  variance  of  teleseismic  magnitude  or  log  amplitude  using 
these  corrections  is  about  25  percent,  similar  to  the  reduction  of  variance  in  teleseismic 
travel  times.  These  results  are  useful  in  predicting  the  structural  resolution  needed  for 
models  of  other  test  sites  to  be  useful  in  making  corrections  for  focusing  and 
defocusing.  The  NTS  results  suggest  that  meaningful  corrections  can  be  made  if  the 
model  resolves  10  to  20  km.  scale  lengths  down  to  100  km.  with  perturbations  of  velocity 
exceeding  4  percent.  Velocity  inversions  that  primarily  resolve  scale  lengths  larger  than 
these  or  that  smooth  over  anomalies  larger  than  2  percent,  (e.g.,  Montfort  and  Evans, 
1982;  Minster  et  al.,  1981),  are  much  less  useful  in  formulating  amplitude  corrections. 

By  analogy  to  the  Taylor  inversion  for  NTS  structure,  the  data  required  to  resolve 
structure  having  these  scale  lengths  would  consist  of  origin  times  and  locations  of  tests 
widely  distributed  over  linear  dimensions  on  the  order  of  100  km.  with  siginificant 
concentrations  of  tests  spaced  less  than  10  km.  apart.  Jt  is  also  necessary  to  obtain  an 
average  crustal  structure  within  the  test  site  from  seismograms  recorded  at  local  and 
regional  range. 

The  focusing  and  defocusing  effects  of  20  to  50  km.  scale  length  structure  having 
perturbation  in  P  velocity  of  several  percent  is  nearly  independent  in  the  frequency 
band  of  teleseismic  body  waves.  This  conclusion  is  even  stronger  in  the  0.2  to  10  Hz. 
band  in  which  measurements  are  made  on  the  teleseismic  P  waves  of  underground 
nuclear  tests.  Frequency  dependent  effects  in  the  coda  of  teleseismic  P  waves  are 
probably  due  to  either  the  effects  of  heterogeneities  having  scale  lengths  smaller  than 
20  km.  and/or  to  frequency  dependent  effects  in  the  scattering  processes  occurring 
near  the  source  and  receiver. 

Magnitude  measurements  based  on  the  integrated  energy  in  the  coda  of  teleseismic 
P  waves  are  likely  to  be  more  stable  because  they  can  remove  some  of  the 
focusing/defocusing  effects  of  three-dimensional  mantle  structure  near  the  source.  The 
deeper  in  the  coda,  the  measurement  is  made,  the  less  affi  led  it  will  be  by  mantle 


structure  near  the  source.  The  optimal  time  in  the  coda  for  this  measurement  should 
be  as  long  as  possible  after  the  direct  P  wave  given  the  signal  to  noise  ratio.  The 
minimum  time  to  achieve  good  stability  can  be  estimated  by  dividing  the  length  of 
characteristic  scale  lengths  of  mantle  structure  near  the  source  by  the  velocity  of  the 
presumed  scattered  wave  near  the  source.  For  example,  assuming  a  3.3  km./sec.  S 
wave  is  scattered  into  a  P  wave  that  is  propagated  to  teleseismic  range,  one  would 
estimate  that  after  30  sec.  into  the  coda,  the  focusing/defocusing  effects  of  100  km. 
scale  and  smaller  length  structure  beneath  the  source  would  be  minimized.  Coda 
magnitudes,  however,  are  only  partially  successful  in  removing  the  focusing/defocusing 
effects  of  structure  beneath  the  source.  They  cannot  remove  these  effects  from  the 
fraction  of  the  coda  that  is  due  to  scattering  of  direct  P  near  the  receiver. 

These  conclusions  are  consistent  with  tests  of  the  relative  performance  of  coda 
versus  classical  magnitudes.  The  predicted  behavior  of  coda  amplitude  critically 
depends  on  assumptions  about  the  distribution  of  mantle  heterogeneity  with  depth. 
Smaller  scale  heterogeneities  with  greater  percent  velocity  fluctuations  are  assumed  to 
be  concentrated  closer  to  the  surface.  Scale  lengths  on  the  order  of  severed  kilometers 
to  10  kilometers  are  assumed  in  the  crust  and  scale  lengths  on  ther  order  of  10  to  100 
kilometers  are  assumed  in  the  upper  mantle.  Fluctuations  in  the  mid  and  lower  mantle 
down  to  the  D"  layer  near  the  core  are  assumed  to  be  smaller  than  1  percent. 

12  2  Descending  Slab  Structures 

A  P  velocity  model  of  the  Kuril-Kamchatka  lithospheric  slab  determined  from  travel 
time  study  by  Creager  and  Jordan  (1986)  has  been  parameterized  in  three  dimensions  to 
investigate  amplitude  and  waveform  effects  on  body  waves.  Very  broad  Gaussian  beams 
were  used  in  a  reciprocal  sense,  shooting  from  receiver  to  source,  to  synthesize  broad 
band  S  waves.  The  per  cent  velocity  perturbation  on  S  velocity  was  assumed  to  be  equal 
to  that  in  P  velocity.  Frequency  dependent  effects  in  the  tail  of  the  S  pulse  were 


observed  in  all  azimuths  on  the  side  of  the  slab  dipping  away  from  the  arc.  This  ‘‘slab 
diffracted"  phase  rapidly  decays  away  within  10  to  15  degrees  of  the  azimuth  parallel  to 
the  strike.  The  character  of  the  slab  diffraction  agrees  with  that  seen  in  finite  difference 
calculations  by  Vidale  (1986)  for  P  waves  propagating  down  dip  from  an  hypocenter  in 
the  slab. 

The  results  of  these  studies  will  be  useful  in  interpreting  the  waveform  studies  of  P 
and  S  attenuation  of  the  type  reported  by  Choy  and  Cormier  (1996),  as  well  as 
investigating  the  shadowing  and  siab  multipathing  evident  in  teleseismic  P  waves  from 
U.S.  nuclear  tests  on  Amchitka  Island  in  the  Aleutian  arc  (Davies  and  Julian,  1972). 

2  Publications 

The  detailed  results  of  the  research  conducted  under  this  contract  are  reported  in 
two  referreed  publications: 

Cormier,  V.F.,  An  application  of  the  propagator  matrix  of  dynamic  ray  tracing:  The 

focusing  and  defocusing  of  body  waves  by  three-dimensional  velocity  structure 
in  the  source  region,  Geophys  J  R  Astr  Soc  ,  87,  1159-1180.  1986. 

Cormier,  V.F.,  Focusing  and  defocusing  of  teleseismic  P  waves  by  known  3-D  structure 

beneath  Pahute  Mesa,  Nevada  Test  Site,  Bull  Seism  Soc  Am  ,  77,  in  press  (will 
appear  in  October  issue),  1987. 

The  first  two  papers  are  included  in  reprint  form  in  sections  3  and  4.  Results  of  the 
slab  modeling  are  still  being  prepared  for  publication.  Preliminary  results  and  modeling 
are  included  in  a  section  5  following  the  reprints  of  the  two  papers  referenced  above. 


3  An  application  of  the  Propagator  matrix  .... 


(Reprint  of  a  paper  by  the  P.I.,  Geophys.  J.  R.  Astr.  Soc.,  87,  1159-1180,  1986.) 
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An  Application  of  the  Propagator  Matrix  of  Dynamic  Ray  Tracing: 
The  Focussing  and  Defocussing  of  Body  Waves  by  Three-Dimensional 
Velocity  Structure  in  the  Source  Region 


V.F.  Cormier 


Earth  Resources  Laboratory 

Department  of  Earth,  Atmospheric,  and  Planetary  Sciences 
Massachusetts  Institute  of  Technology 
Cambridge.  MA  02139 


Geophys  J  R  Asfr  Soc  .  87.  1936.  1169-1190 


Summary.  Since  the  dynamic  ray  tracing  system  can  be  written  in  the 
linear  form  =  SW  ,  it  permits  the  definition  of  fundamental  matrix  and 

propagator  matrix  solutions.  The  propagator  matrix  can  be  exploited  to 
connect  complex  3-D  and  2-D  regions  to  1-D  regions  of  a  modeL  For  example,  if 
the  three-dimensional  variations  of  a  model  are  confined  only  to  the  portion  of 
the  ray  path  between  0  and  00,  numerical  integration  of  the  linear  system  need 
only  be  performed  between  0  and  0a  to  construct  the  propagator  II (0a,0)  . 
Along  the  segment  of  the  ray  path  between  O'  and  C^,  traversing  a  1-D  or 
homogeneous  portion  of  the  model,  analytic  solutions  exist  for  11(0, Q)  .  Thus 
n(00,Og)  can  be  simply  obtained  by  multiplying  the  analytic  and  numerically 
obtained  II  matrices.  With  this  procedure,  quantities  are  calculated  that  are 
needed  for  the  synthesis  of  teleseismic  P  waves  by  superposition  of  Gaussian 
beams.  These  synthetics  are  used  to  investigate  the  focussing  and  defocussing 
of  teleseismic  P  waves  by  3-D  structure  in  the  vicinity  of  the  source.  A  source 
region,  having  a  fluctuation  in  velocity  of  4  per  cent  over  a  characteristic  scale 
length  of  50  to  100  km.,  produces  factors  of  two  fluctuations  in  amplitude  and 
several  tenths  of  a  second  in  travel  time  at  teleseismic  range  in  experiments  in 
which  either  source  location  is  varied  at  constant  azimuth  or  in  which  azimuth 
is  varied  at  constant  source  location.  A  model  having  a  maximum  fluctuation  as 
small  as  0.8  per  cent  is  capable  of  producing  caustics  and  multipaths  at 
teleseismic  range,  depending  on  its  distribution  of  scale  lengths  and  its  ratio  of 
vertical  with  respect  to  horizontal  scale  length.  The  multipaths  and  caustics  of 
such  a  model,  however,  occur  over  too  small  an  area  and  are  too  closely  spaced 
in  arrival  time  to  be  resolved  with  seismograph  systems  in  the  0.01  to  4  Hz. 


1  Introduction 


All  of  the  asymptotically  approximate  methods  of  synthesizing  body  waves 
in  three-dimensionaily  varying  media  are  closely  related.  They  all  require  the 
calculation  of  travel  times  and  either  geometric  spreading  or  quantities  closely 
related  to  geometric  spreading.  In  the  case  of  ray  theory,  a  travel  time  and 
spreading  factor  must  be  estimated  between  two  given  points.  In  the  case  of 
Gaussian  beams  (terven^,  1983)  or  WKBJ/Maslov  seismograms  (Chapman  and 
Drummond,  1982)  travel  times  and  weighting  factors  that  depend  on  ray 
spreading  must  be  calculated.  The  weighting  factors  appear  in  the  integrand  of 
a  superposition  integral  over  either  take-off  angles  or  slowness  components 
measured  at  the  source  point.  In  applications  of  the  Kirchhoff-Helmholtz 
integral  (e.g.,  Haddon  and  Buchen,  1981;  Scott  and  Helmberger.  1983),  travel 
time  and  geometric  spreading  must  be  estimated  at  points  on  a  spatial 
integration  surface. 

In  any  of  these  techniques,  the  travel  times,  weighting  factors,  and 
geometric  spreading  can  be  calculated  in  a  three-dimensionally  varying  model 
by  numerical  integration  of  linear  differential  equations.  These  equations 
consist  of  a  kinematic  system  for  ray  trajectory  and  travel  time  (e.g.,  terven^ 
et  al.,  1977)  and  a  dynamic  system  (Cerven^  and  Hron,  1980)  for  quantities 
related  to  the  weighting  factors  and  spreading,  The  kinematic  system  requires 
the  specification  of  seismic  velocity  and  its  first  spatial  derivative.  In  addition 
to  these,  the  dynamic  system  requires  the  specification  of  the  second  order 
spatial  derivatives  of  velocity. 

The  computational  expense  of  numerical  integration  of  the  kinematic  and 
dynamic  systems  increases  with  the  length  of  ray  paths.  Hence,  it  is  desirable 
to  incorporate  the  known  analytic  forms  of  the  integrated  equations  along 
portions  of  the  ray  path  traversing  regions  of  the  model  that  are  either 


homogeneous  or  vary  in  only  one  co-ordinate  direction.  The  patching  of 
analytic  to  numerical  solutions  is  particularly  useful  in  the  synthesis  of  body 
waves  at  teleseismic  distances.  In  this  case,  long  segments  of  the  ray  paths  are 
contained  in  the  mid-mantle  between  700  km.  and  2500  km.  depth,  where  all 
evidence  suggests  that  velocity  fluctuations  and  departures  from  radial 
symmetry  are  much  weaker  than  either  in  the  crust  and  upper  mantle  or  in  the 
lowermost  mantle  (e.g.,  Dziewonski,  1984;  Woodhouse  and  Dziewonski,  1984). 

In  the  following  sections,  a  procedure  is  described  for  the  multiplying  of 
propagator  matrices  to  patch  analtyic  to  numerical  solutions  of  the  equations 
of  dynamic  ray  tracing.  The  analytic  forms  for  the  elements  of  this  propagator 
are  given  for  a  flattened.  1-D  Earth  model.  Calculations  using  this  propagator 
are  illustrated  for  the  synthesis  of  P  waves  at  teleseismic  distances  by 
superposition  of  Gaussian  beams.  The  specific  problem  considered  in  these 
calculations  is  the  focussing  and  defocussing  of  body  waves  due  to  3-D 
structure  in  the  region  of  the  source.  One  example  of  such  focussing  and 
defocussing  is  the  shadowing  effect  of  a  subducted  lithospheric  slab  (Sleep, 
1973).  Other  examples,  which  are  discussed  in  this  paper,  are  related  to  yield 
estimation  of  underground  nuclear  explosions.  At  a  constant  distance,  the 
focussing  and  defocussing  effects  of  3-D  structure  in  the  source  region  can 
introduce  an  azimuthal  variation  in  the  amplitudes  of  P  waves  radiated  from  an 
explosive,  isotropic  source.  This  focussing  and  defocussing  also  introduces 
variations  in  the  apparent  yield  of  explosions  of  equal  yield  as  their  location 
varies  within  a  test  site. 

2  The  linear  equations  of  dynamic  ray- tracing 

Both  travel  times  and  dynamic  spreading  quantities  can  be  estimated  in  3- 
D  elastic  media  from  a  paraxial  approximation  of  the  wavefront,  in  which  2x2 


matrices  P  and  Q  are  defined  (  Cerven^  ,  and  Psencik  ,  1983  )  .  At  any  given 
point  along  a  ray  the  values  of  P  and  Q  can  be  found  by  integrating  a  system  of 
linear  equations: 

^■  =  kP  (la) 

as 

and 


^-  =  ^-*VQ  (lb) 

Here  ds  is  an  incremental  element  of  ray  path,  v  is  a  local  velocity,  and  V  is  a 
matrix  of  second  derivatives  with  respect  to  ray  centered  co-ordinates  that 
move  along  the  ray.  The  elements  of  V  are 


i/u  u,2 

^  [u21  i>22 

The  subscripts  on  the  elements  of  V  denote  differentiation  of  velocity  v  with 
respect  to  the  ray  centered  co-ordinate  directions  qx  and  q2 

Equations  (la-b)  can  be  rewritten  as  a  single  linear  system  consisting  of 
four  equations: 


(2) 


£  =  SW  (3) 

where  Vis  a  column  vector  having  four  components  and  S  is  a  4  x  4  matrix.  The 
form  of  such  a  system  is  similar  to  the  system  for  stress-displacement 
components  in  a  vertically  varying  medium.  Here,  of  course,  the  differentiation 
is  with  respect  to  incremental  path  length  along  a  ray  rather  than  depth.  Just 
as  for  the  stress-displacement  system,  propagator  matrix  solutions  can  be 
found  to  this  system.  Specifically,  Cerven^  (1985a)  has  defined  a  fundamental 
matrix  for  the  equation  system  (3).  which  can  be  normalized  to  a  4x4 
propagator  matrix  T[(00,0, )  by  choice  of  initial  conditions  at  a  point  0t  ,  so  that 
W  (00)  =  WO,  0t  )W(0f )  .  Since  FI  is  a  propagator  matrix,  it  has  the  property  that 
IX  o.  ,0,)  =  i  .  where  I  is  the  identity  matrix.  Also,  given  a  point  0  between  Q, 
and  0„  .  fl (0g,0,)  =  W{00,O)  IT (0.0,)  .  This  property  can  be  used  to  accelerate 


two  point  ray  tracing  and  the  calculation  of  spreading  functions  and  weighting 
factors  needed  for  synthesis  of  body  waves  in  three-dimensionally  varying 
media.  For  example,  if  the  three-dimensional  variations  of  a  model  are  confined 
only  to  the  portion  of  the  ray  path  between  0  and  09,  numerical  integration  of 
the  linear  system  need  only  be  performed  between  0  and  09  to  construct 
r \09,0)  .  Along  the  segment  of  the  ray  path  between  0  and  Og,  traversing  a  1- 
D  or  homogeneous  portion  of  the  model,  analytic  solutions  exist  for  n (0,0,) 
(Figure  1).  Thus  I1(09,H)  can  be  simply  obtained  by  multiplying  the  analytic 
and  numerically  obtained  II  matrices. 


3  Hie  elements  of  the  propagator  matrix 


3.1  GENERAL  3-D 


The  II  matrix  can  be  defined  in  terms  of  2  x  2  sub  matrices  P  and  Q  as 
follows: 


n  = 


Q 1  Q* 

pi  pff 


(4) 


The  superscripts  I  and  R  distinguish  solutions  to  the  dynamic  ray  tracing 
equations  (la-b)  using  plane  wave  and  point  source  initial  conditions 
respectively.  The  initial  conditions  for  Q^,  P^  and  Q*.  P*  have  been  discussed 
extensively  by  terven^  and  Psencik  (1983a, b)  and  Cerven^  (1985a, b)  in  the 
context  of  the  general  complex  solutions  of  eq.  (3)  for  Gaussian  beams.  The 
simplest  form  of  the  initial  conditions  are  that  at  0, 

Q*  =  0  ;  P*  =  I 


(5a) 


and 


O'  =  I  ;  P'  =  0  (5b) 

Eqs.  (5a, b)  are  necessary  for  II  to  have  propagator  character;  more  general 

conditions  can  be  incorporated  by  specifying  W(0f )  . 


3.2  1-D 


& 


In  order  to  determine  the  elements  of  the  TT  matrix,  a  vector  basis  must 
first  be  selected  for  the  ray  centered  co-ordinate  system  (  gi.?2.s  )  at  the  point 
at  which  initial  conditions  are  specified.  Here,  given  a  fixed  Cartesian  co¬ 
ordinate  system (x.  y,  z),  the  ray  will  be  assumed  to  be  confined  to  the  (x,  z) 
plane  in  a  1-D  medium  and  the  e]  direction  will  be  chosen  to  point  in  the  y 
direction.  A  vector  eg  is  chosen  such  that  (  ei  eg,  t)  forms  aun  orthogonal, 
right-handed  vector  basis,  where  t  is  the  ray  tangent. 

With  this  definition  of  the  ray  centered  co-ordinate  system,  and  assuming 
that  velocity  varies  only  as  a  function  of  z  ( 1-D  medium),  it  is  possible  to  derive 
analytic  expressions  for  the  elements  of  R  The  solutions  for  the  spreading  of  a 
line  source  in  a  vertically  (z)  varying  medium  have  been  given  by  Madariaga 
(1984).  These  are  appropriate  for  the  components  Qy  .  Py  .  i.j  =  2,2  .  For  the 
spreading  of  a  point  source  in  three  dimensions,  solutions  for  the  components 
i.j  =1,1  must  be  added  to  these.  The  components  i.j  =1.1  are  quite  simple  in 
form.  This  is  because  by  the  choice  of  the  ray  centered  co-ordinates  in  the  1-D 
medium  =0.  The  off  diagonal  elements  of  the  Q  and  P  matrices  are  all  zero 
because  the  velocity  v  is  assumed  not  vary  in  the  x  and  y  directions.  In 
summary,  the  16  components  of  the  FI  matrix  in  a  vertically  varying  medium  are 
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If  a  ray  has  a  turning  point,  then  the  angle  6  must  be  in  the  range  0  <  6  <n  . 
When  a  ray  is  propagating  downward  in  the  direction  of  the  positive  depth  (z) 

axis,  <5  is  in  the  range  (0,  ^-)  ;  when  it  is  propagating  upward,  6  is  in  the  range 

(0,rr)  .  This  extended  angular  range  is  important  for  properly  describing  plane 
wave  caustics  at  turning  points  (Madariaga,  1984).  It  is  also  essential  to  include 
in  the  calculation  of  IT  elements  in  a  product  of  II  matrices  whenever  one  of  the 
member  matrices  is  associated  with  a  turning  ray. 


3.3  CONNECTING  REGIONS 


With  the  elements  of  the  propagator  n  above,  it  can  be  verified  that  for  any 
point  O  along  the  ray  path  at  which  the  velocity  and  its  first  spatial  derivatives 
are  continuous,  that 


x(oa)  =  n(  o0.  o  )ii(  o'.  q,)-x(o,)  (7) 

where  X(0a)  is  a  4  x  2  matrix  of  initial  conditions  on  P  and  Q ,  and  X(0„)  is  a  4  x 
2  matrix  formed  from  a  linear  combination  of  plane  wave  (1)  and  point  source 


solutions  (R)  for  P  and  Q  .  For  example,  point  source  initial  conditions  at  0, 


can  be  specified  by  setting  X(C^) 


0  0 
0  0 
1  0 
0  1 


Then  by  eqs.  (7)  and  (4).  X(0„) 


Qff 

P* 


It  can  also  be  shown  by  substitution  of  the  definitions  of  eqs.  (6),  that 
IXO'.Q.)  =  n (O0,O)  IT (0,0,)  ,  that  the  determinant  of  II  is  constant  along  the 
ray  path,  and  that  det  n  =  1 . 

If  the  velocity  or  its  first  spatial  derivative  are  discontinous  at  0  ,  then  the 
jump  conditions  defined  by  CervenJ’  (1985a)  must  first  be  applied  to  the 
elements  of  ri(0,0,)  before  the  multiplication  of  eq.  (7)  is  performed.  These 
jump  conditions  are  identical  for  both  plane  wave  and  point  source  solutions. 
The  jump  conditions  can  be  compactly  expressed  by  a  4  x  4  matrix  F,  with  eq. 

(7)  becoming 

no0)  =  U(o0.o*)rn(o-,a.)x(o.)  (0) 

where  the  +  sign  refers  to  quantities  evaluated  on  the  incident  side  of  the 
boundary  at  O  and  the  -  sign  refers  to  quantities  evaluated  on  the  transmitted 
side  of  the  boundary  at  O'.  The  elements  of  F  are  given  in  Cerven^  (1985a)  for 
the  most  general  case  of  discontinuities  in  velocity  and  its  first  spatial 
derivatives.  For  the  co-ordinate  system  chosen  here  and  for  velocity  depending 
only  on  z,  the  elements  of  F  are  all  zero  except  for 


O) 


Thus,  computations  with  the  propagator  matrix  of  dynamic  ray  tracing  are  not 


completely  analagous  to  those  using  the  stress-displacement  propagator.  In 
the  former  case,  the  dynamic  ray  tracing  propagator  is  discontinuous  at  first 
and  second  order  velocity  discontinuities,  while  in  the  latter  case  the  stress- 
displacement  propagator  is  continuous  at  solid- solid  boundaries.  The  patching 
of  solutions  in  1-D  and  homogeneous  regions  to  those  in  3-D  regions  is  most 
conveniently  performed  at  pseudo  boundaries  at  which  velocity  and  its  first 
spatial  derivatives  are  continuous.  Patching  can  be  performed  at  more  general 
boundaries,  but  the  jump  conditions  must  first  be  used  to  correct  P  and  Q 
elements  to  values  appropriate  for  the  transmitted  side  of  the  boundary. 

A  final  note  should  be  made  that  in  order  for  the  II  matrix  to  have  the 
properties  of  a  propagator,  one  must  use  the  plane  wave  initial  conditions  (5b) 
in  calculating  P^,  rather  than  the  modified  plane  wave  conditions 
recommended  by  Madariaga  (1994)  for  a  source  in  a  region  having  a  non-zero 
velocity  gradient.  The  modified  initial  conditions  of  Madariaga  (1984)  were 
proposed  to  obtain  phase  fronts  at  the  source  that  were  equivalent  to  plane 
waves.  This  condition  was  necessary  to  stabilize  the  computation  of  Gaussian 
beam  seismograms.  It  also  makes  Gaussian  beam  synthesis  more  closely 
resemble  the  process  of  plane  wave  superposition  upon  which  the  WKBJ/Maslov 
technique  (Chapman  and  Drummond,  1982)  is  formulated.  Recently,  however, 
Cerveny'  ( 1985a, b)  has  proposed  optimal  beam  width  parameters  that  are 
equivalent  to  plane  wave  fronts  at  the  receiver.  These  parameters  are  optimal 
in  the  sense  of  minimizing  the  error  in  discretizmg  and  truncating  the  integral 
that  superposes  beams.  This  choice  of  beam  parameters  also  addresses  the 
stability  problem  noted  by  Madariaga.  For  a  source  and  receiver  both  at  the 
surface,  terveny"s  optimal  beam  width  parameters  are  reciprocally  equivalent 
to  Madariaga's  modified  plane  wave  initial  conditions.  This  can  be  demonstrated 
using  the  propagator  matrix  of  dynamic  ray  tracing  and  reciprocal  relations 
between  elements  of  the-P  and  Q  matrices. 


Section  4,  which  follows,  reviews  the  procedures  used  in  synthesizing  the 
seismograms  shown  and  discussed  in  section  5  .  These  P  wave  seismograms 
were  synthesized  by  a  summation  of  Gaussian  beams,  in  which  the  propagator  FI 
was  used  to  connect  a  region  having  3-D  variations  in  velocity  to  a  1-D  region. 

4  Gaussian  beam  summation  in  a  3-D  model 

4.1  THE  TECHNIQUE 

Seismogram  synthesis  by  Gaussian  beams  in  three-dimensionally  varying 
media  has  been  fully  described  by  Cerven^  ( 1985a, b).  The  steps  are  repeated 
here  only  in  summary  form,  emphasizing  the  particular  choices  made  for  the 
parameters  affecting  beam  widths  and  special  considerations  in  the  use  of 
propagator  matrices. 

A  spectral  component  of  vertical  displacement  observed  at  point  5  was 
calculated  by 

u(u,  S)  =  J f  *N(yx.72)0N(O,)exp[i'MS.  00)  -  l^-k(0,.00)]dy^  (10) 

Seismograms  were  synthesized  in  the  time  domain  by  the  spectral  method,  in 
which  eq.  (10)  was  first  evaluated  at  the  discrete  frequencies  required  by  a  fast 
Fourier  transform.  Source  spectra  and  seismograph  responses  were  applied  in 
the  frequency  domain  before  inverse  transforming  to  the  time  domain. 

In  eq.  (10),  ?),  y2  are  take-off  angles  at  the  source  point  0,  .  Vn{7\<yz)  is  a 
weighting  function  that  varies  for  each  set  of  take-off  angles  .  tf*  is  a  a  factor 
that  includes  radiation  pattern  and  source  normalization  but  does  not  include 
geometric  spreading  00  is  the  orthogonal  projection  of  point  S  onto  the  ray. 
For  the  applications  here,  00  can  be  taken  as  the  ray  end  point  on  the  surface 
of  the  Earth  r (5.  )  is  the  Gaussian  beam  complex-valued  travel  time  at  S, 

which  is  defined  m  terms  of  ray  centered  co-ordinates  by 


(11) 


t(S,0o)  =r(09)  +  |-qrM(00)q 
M  is  a  2  x  2  symmetric  matrix  of  the  second  derivatives  of  the  complex  valued 
travel  time  with  respect  to  the  ray  centered  co-ordinates  qr  =  (g,f  gz).  U  is 
evaluated  from 

U  =  PQ~l  (12) 

k(0t,0a)  is  an  integer  index,  which  Ziolkowski  and  Deschamps  (1980)  have 

termed  the  KMAH  index.  The  KMAH  index  specifies  the  number  of  phase 

advances  encountered  as  each  caustic  is  passed  by  a  ray.  Practical  calculation 
of  the  KMAH  in  three-dimensionally  varying  media  is  discussed  in  greater  detail 
in  Chs.  4  and  6  of  Cerven^  ( 1985a).  Examples  of  its  calculation  in  2-D  media  are 
given  in  Chapman  and  Drummond  (1982),  Nowack  and  Aki  (1984),  and  Cormier 
and  Spudich  (1984). 


i 
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4.2  BEAM  PARAMETERS 


The  P  and  Q  matrices  of  eq.  (12)  are  generally  complex.  The  complex 
matrices  P  and  Q  are  completely  specified  at  any  point  by  the  2x4  matrix  X  of 
eq.  (7),  with  P  and  Q  being  the  lower  and  upper  2x2  submatrices  respectively 
of  the  X  matrix.  Cerven^  ( 1985a, b)  specifies  beam  parameters  using  M  rather 
than  P  and  Q  individually.  This  is  analagous  to  imposing  initial  conditions  on 
the  "impedance  ratio"  X  =  XQ'at  ray  end  points.  The  propagator  relation  of  eq. 
(7)  makes  it  possible  to  specify  complex  initial  conditions  at  any  point  along  a 
ray  path.  Thus  at  the  ray  end  point  Og,  X (00)  can  be  specified  by 

*(q,)=p<J-.  =  I  (13) 

At  any  other  point  along  the  ray,  say  point  0,  at  the  source, 

X(Ot)  =  n(0',00)X(00)  (14) 

In  the  example  calulations,  the  initial  conditions  were  specified  on  ii  at  the  ray 


I $3 


end  points,  as  in  eq.  (13).  The  source  point  0,  and  ray  end  points  00  were 
chosen  to  be  in  locally  homogeneous  layers.  Re  M (00)  was  specified  by 
Cerven^'s  "effective  plane-wave"  condition.  For  a  homogeneous  layer,  this 
condition  reduces  to 


ReU  =  0  (15) 

where  0  is  a  2  x  2  matrix  of  zeros.  The  effective  plane  wave  condition  is  the  set 

of  initial  conditions  that  produce  zero  second  derivatives  of  the  travel  time  field 

along  the  Earth’s  surface.  With  this  choice,  the  real  part  of  the  complex  travel 

time  t(S,  00)  becomes 

r(S,00)  =T(0a)+p-[(*(S)-*(09)]  (16) 

where  p  is  the  vector  slowness  at  the  end  point  of  the  ray;  x(S)  is  the  cartesian 

position  vector  of  the  receiver;  and  x(00)  is  the  cartesian  position  vector  of  the 

ray  end  point.  In  this  form,  the  real  part  of  the  complex  travel  time  t(S,00) 

resembles  the  phase  function  in  the  WKBJ/Maslov  technique  (Chapman  and 

Drummond,  1982). 

Im  H  was  selected  at  ray  end  points  to  minimize  the  error  in  the 
discretiztion  of  the  superposition  integral.  If  the  ray  end  point  is  in  a 
homogeneous  layer,  this  condition  reduces  to 

Im  II  =  iH*|.  (17) 

where  denotes  an  absolute  value  of  the  matrix  II*  =  P*(Q*)  1  •11*1  can  be 

calculated  by  finding  the  matrix  of  0  of  eigenvectors  of  II*  associated  with 

eigenvalues  A^Ag  and  computing 


i  ii*  =  n 


^i  I 

o 


o 

xz  | 


n7 


(18) 


Note  that  the  choices  ( 15)  and  ( 17)  for  complex  M  at  the  ray  end  point  00 
uniquely  determine  the  value  of  complex  M  at  the  source  or  starting  point  of 
integration  0 ,  The  value  of  complex  II  at  0 ,  may  be  found  by  propagator 
multiplication  as  in  eq  (14)  The  value  11(0,)  ,  however,  does  not  need  to  be 
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calculated  in  order  to  evaluate  the  integrand  of  eq.  (10).  Moreover,  beam 
parameters  may  be  specified  at  any  point  along  a  ray  through  a  complex  li 
matrix.  For  example,  one  may  specify  the  parameters  at  the  point  where  a  ray 
is  incident  on  some  boundary.  The  values  of  complex  M,  which  control  beam 
decay  and  interference,  are  then  uniquely  determined  at  either  end  point  of 
the  ray  by  the  propagator  relation. 

If  initial  conditions  are  specified  at  the  end  points  00l  the  weighting 
function  tfiN  in  eq.  (9)  becomes 

**  =  ^  !  de  tQ^  |  [  -det(M(  0o )  -M*(  0„ )  ]  ^  (19) 

with  the  branch  cut  of  the  square  root  taken  such  that 

1 

Re[  -det(lf(  0„ )  — M* ( ) )  ] 2  >0. 

The  amplitudes  U*  were  taken  to  be  constant  over  all  beams  and  source 
locations  in  the  example  calculations  described  in  the  next  section.  This 
neglects  effects  of  the  additional  component  of  Gaussian  beams,  as  well  as 
factors  having  to  do  with  source  normalization  and  the  free  surface  coefficient. 


In  all  examples  shown,  both  the  additional  component  and  the  free  surface 
coefficient  had  insignificant  variations  over  the  spot  ellipse  (  the  boundary 
surrounding  a  beam  at  which  e  decay  is  achieved).  Likewise  in  experiments 
involving  variations  in  source  locations,  velocity  variation  was  too  small  to 
significantly  affect  a  source  normalization  factor. 

With  beam  parameters  specified  at  ray  end  points,  the  evaluation  of  eq. 
(10)  only  requires  knowledge  of  P*,  Q*  at  ray  end  points.  P1 ,  ,  however,  must 

also  be  known  along  both  analytically  and  numerically  integrated  segments  if 
propagator  products  are  evaluated. 

The  calculation  of  the  KMAH  index  is  simplified  because  any  square  roots 
appearing  in  the  integrand  of  eq.  (10)  depend  only  on  P*.  Q^.  The  changes  in 
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sign  of  det  Qff 

and  the  trace  of  Q *  are  tracked  along  rays  in  the  numerically 
integrated  regions  of  a  model.  Along  the  analytically  integrated  sections  of  a 
ray  path,  the  signs  of  these  quantities  are  compared  before  and  after  the 
propagator  multiplications.  The  KMAH  index  is  advanced  one  unit  if  the  sign  of 
det  Q"  changes  but  the  sign  of  the  trace  of  remains  unchanged.  The  KMAH 
index  is  advanced  two  units  if  the  signs  of  both  det  Q*  and  the  trace  of 
change  (  see  Cerven^,  1985a)  .  For  this  procedure  to  work  along  segments  that 
are  analtyically  integrated  by  propagator  multiplication,  it  is  important  to  allow 
the  angle  <5  ,  which  appears  in  all  analytically  computed  elements,  to  be  defined 

over  an  angular  range  that  exceeds  whenever  a  ray  passes  through  a  turning 
point. 

4.3  BEAM  DECAY  AND  INTERFERENCE 

Figures  2  and  3  illustrate  how  the  waveforms  obtained  from  Gaussian  beam 
superposition  critically  depend  on  two  phenomena:  the  exponential  decay  rate 
away  from  the  central  ray.  and  the  constructive  interference  of  the  phase  of 
neighboring  beams.  In  Figure  2,  the  values  of  the  exponential  decay  power  at  1 
Hz.  are  plotted  at  the  beam  end  points  in  the  vicinity  of  a  teleseismic  station  at 
70°  distance  for  beams  departing  from  one  of  the  3-D  models  of  a  source  region 
that  are  discussed  in  the  next  section.  These  decay  rates  were  calculated  for 
the  beam  parameters  given  by  eqs.  ( 15)  and  (17).  Note  that  a  large  region  can 
potentially  contribute  to  a  1-Hz.  signal.  Lower  frequencies  would  have  a 
proportionally  smaller  beam  decay.  A  0.03  Hz  signal  would  have  contributions 
from  beams  as  much  as  1000  km  or  more  away  from  the  station.  Equally 
important,  however,  are  the  effects  of  constructive  interference  shown  in 
Figure  3.  For  these  beam  parameters,  the  region  of  constructive  interference 
of  a  1  Hz  wavelet  is  slightly  larger  than  the  region  in  which  exponential  decay  is 
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5  Results  with  a  3-D  model  of  the  source  region 

5. 1  A  MODEL  INVERTED  FROM  TRAVEL  TIMES 

The  focussing  and  defocussing  of  teleseismic  body  waves  by  3-D  structure 
in  the  vicinity  of  the  source  have  been  investigated  with  two  different  models. 
The  first  model  is  one  obtained  by  Zandt  (1981)  for  central  California  using  a 
block  inversion  of  teleseismic  travel  times  by  the  method  of  Aki  et  al.  (1976). 

The  Zandt  (1981)  model  has  four  layers  from  0.0  to  90.0  km  in  depth.  The 
horizontal  block  size  is  10.0  km  in  the  top  layer  and  20.0  to  25.0  km  in  the  lower 
layers.  Average  velocity  variations  are  between  4.0  to  8.0  percent  in  the  top 
layer  and  2.0  to  4.0  percent  in  the  lower  layers.  The  rms  velocity  variation 
measured  over  grid  points,  however,  is  generally  much  lower,  on  the  order  of 
less  than  2  percent.  This  is  because  the  largest  variations  take  place  over 
relatively  broad  regions,  having  characteristic  scale  lengths  of  between  50  to 
100  km.  (Figure  4). 

Seismograms  were  synthesized  in  this  model  by  summation  of  Gaussian 
beams.  An  explosive  point  source  was  assumed  at  the  center  of  the  model  at  9.6 
km.  depth.  The  FI  elements  were  calculated  in  the  3-D  source  region  by 
numerical  integration  of  the  kinematic  and  dynamic  ray  tracing  equations. 
Velocities  and  their  first  and  second  order  spatial  derivatives  in  the  3-D  region 
were  defined  by  the  coefficients  of  cubic  spline  interpolators  between  grid 
points.  The  3-D  region  was  patched  into  a  1-D,  radially  symmetric  Earth  model 
at  90  km.  depth.  The  1-D  Earth  model  was  the  1  Hz.,  isotropic  PREM  of 
Dziewonski  and  Anderson  (1981).  PREM  was  first  flattened  using  the 
transformations  described  by  Muller  (1971)  .  The  TI  elements  of  eq.  (6)  were 
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Figure  4.  P  velocity  contours  in  a  horizontal  plane  at  90  km.  depth  in  the  3-D 
model  for  central  California  by  Zandt  ( 1991).  Also  shown  are  the  projec¬ 
tions  of  source  locations  used  in  propagation  experiments  at  teleseismic 


then  computed  using  a  fast  algorithm  in  which  the  quantities  — —  and  X  are 

given  by  a  analytic  formulae  summed  over  thick,  vertically  inhomogeneous 
layers  (Cerven^  and  Jansky.  1983).  n  was  determined  at  the  ray  end  points  in 
the  receiver  region  by  propagator  multiplication.  The  matrices  P^,  Q*  needed 
to  evaluate  eq.  (10)  are  given  by  the  2x2  sub-matrices  of  FI ,  rta  and  ni2 
respectively.  Focussing/defocussing  effects  of  the  structure  were  calculated  at 
teleseismic  range  for  the  variations  in  azimuth  and  variations  in  lateral  source 
location  shown  in  Figure  4. 

5. 1  1  Effect  of  varying  azimuth,  at  constant  source  location 

Figure  5  shows  the  results  of  beam  summation  for  a  teleseismic  P  wave 
from  an  explosive  source  embedded  in  the  Zandt  (1981)  model.  At  source 
location  SO,  seismograms  were  computed  at  70°  for  eight  different  azimuths.  In 
each  column  of  Figure  5,  the  amplitudes  predicted  in  different  pass  bands  are 
shown.  The  broadband  pulse  was  that  obtained  using  the  source-time  function 
of  Madariaga  and  Papadimitriou  (1985)  in  a  frequency  band  between  0.03  to  4 
Hz.  Amplitudes  are  scaled  to  the  maximum  peak  to  peak  amplitude  observed  at 
the  101°  azimuth.  Amplitude  variations  are  of  the  order  of  two  and  travel  time 
variations  on  the  order  of  several  tenths  of  a  second.  The  travel  time  variations 
are  consistent  with  the  focussing/defocussing  effects  --  large  amplitudes 
correlate  with  slow  travel  times  and  small  amplitudes  correlate  with  fast  travel 
times 

Note  that  even  in  an  expanded  time  scale,  it  would  be  difficult  to  accurately 
measure  the  small  travel  time  fluctuations  associated  with  this  model.  In  Figure 
5,  the  negative  peak  of  the  short  period  record  having  an  amplitude  of  100  units 
lies  just  to  the  right  of  the  reference  line  in  the  center  column.  In  the  record 
having  45  units  amplitude,  this  negative  peak  coincides  or  lies  just  to  the  left  of 


rnograms  constructed  by  superposition  of  Gaussian 


the  reference  line.  These  small  differences  in  arrival  time  do  not  bode  well  for 
any  attempt  to  include  amplitude  variations  with  travel  time  variations  in  an 
inversion  for  velocity  structure  using  real  data. 

The  largest  amplitudes  correlate  with  azimuths  in  which  the  beams  sample 
a  low  velocity  anomaly  to  the  southeast  of  source  SO.  This  anomaly  is  a 
significant  feature  in  both  layers  3  and  4  of  the  Zandt  model,  persisting  over  50 
km.  of  depth  in  the  model.  Zandt  interprets  this  feature  as  a  NW-SE  trend  of 
lithospheric  thinning  associated  with  a  fault  zone  that  includes  the  Calaveras. 
Rogers  Creek.  Maacama,  and  Lake  Mountain  faults. 

5  12  Effect  of  varying  source  location  at  constant  azimuth 

For  a  fixed  azimuth,  and  variations  in  source  site  from  positions  S30-  to 
S30+  (Figure  6).  amplitude  variations  are  small.  This  reflects  the  smaller 
differences  in  structure  between  the  regions  sampled  by  the  beams  compared 
to  those  in  the  azimuthal  experiment.  The  velocity  anomalies  in  the  deeper 
layers  are  broad  features  having  scale  lengths  of  50  km.  or  more.  The 
anomalies  in  the  shallower,  crustal  layers  have  smaller  scale  lengths,  but  the 
crustal  layers  are  thin  compared  to  the  total  thickness  of  the  model  and  the 
beams  spend  much  longer  time  in  the  thick  layers  3  and  4.  Thus  the  broad  scale 
lengths  of  the  anomalies  in  layers  3  and  4  have  the  greatest  influence  on 
amplitudes.  This  is  consistent  with  the  large  variations  in  amplitude  shown  in 
the  azimuthal  experiment  (Figure  5)  as  well  as  with  the  smaller  variations  in 
amplitude  due  to  changes  in  receiver  location  over  a  line  having  a  length 
roughly  equal  to  the  scale  length  of  the  broad  anomalies  (Figure  6) 


5  1.3  Frequency  independence  and  its  implications  for  treaty  verification 

The  relative  amplitudes  in  Figures  5  and  6  are  nearly  independent  of  the 
frequency  band  of  observation.  This  result  has  important  consequences  for  the 
yield  estimation  of  underground  nuclear  explosions  by  measurements  of 
classical  body  wave  magnitudes,  m^,  versus  broader  band  measurements  of 
radiated  energy  in  the  time  and  frequency  domain.  If  deep  seated,  broad  scale 
length  (50km  and  greater),  velocity  anomalies  of  2%  or  more  are  a  common 
occurence  in  the  upper  mantle  of  the  Earth,  they  will  act  to  focus  and  defocus 
body  waves  over  a  broad  frequency  band.  The  focussing  and  defocussing 
caused  by  these  broad  anomalies  will  be  independent  of  frequency  and  will  thus 
introduce  a  scatter  in  broader  band  measures  of  radiated  energy  which  will  be 
equivalent  to  that  seen  in  the  narrow  band  m^  measurement.  Focussing  and 
defocussing  by  structure  in  the  source  region  will  also  affect  the  coda  of  P 
waves  if  a  portion  of  this  coda  is  generated  in  the  receiver  region.  These  effects 
may  help  explain  why  broader  band  and  integrated  coda  measures  of  body  wave 
energy  often  do  not  exhibit  any  less  scatter  than  classical  mb  measurements. 

I  he  broadband  and  coda  magnitudes  that  exhibit  the  least  scatter  typically 
have  0.15  to  0.2  standard  deviation  in  units  of  logarithm  of  energy  flux  rate  over 
source  or  receiver  arrays  having  apertures  of  200  km  (Bullitt  and  Cormier. 

1984)  This  corresponds  to  a  factor  of  1.5  to  2  variation  in  the  amplitudes  of 
particle  velocity,  similar  to  that  seen  in  the  synthetic  seismograms  of  Figures  5 
and  6  These  results  suggest  that  knowledge  of  the  broader  scale  length 
velocity  anomalies  beneath  source  and  receiver  sites  may  be  useful  in 
correcting  and  reducing  the  scatter  m  magnitude  estimates  and  hence  the 
uncertainty  in  yield  estimates  of  nuclear  tests. 

The  frequency  independence  of  the  amplitudes  calculated  in  the  Zandt 
model  is  a  characteristic  of  ray-theoretically  predicted  amplitudes  It  suggests 
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that  the  Gaussian  beam  synthesis  should  not  be  necessary  in  order  to 
accurately  calculate  amplitude  variations  due  to  this  receiver  structure.  There 
are  no  caustics  in  this  example  and  there  is  no  particular  advantage  to  using 
Gaussian  beam  superposition  over  asymptotic  ray  theory,  other  than  exploiting 
the  paraxial  estimation  of  travel  time.  Additional  evidence  of  why  beam 
superposition  in  this  example  reproduces  simple  ray  theory  is  illustrated  in 
Figures  7.  Ray  densities  (Figure  7)  within  the  vicinity  of  a  receiver  can 
accurately  predict  the  amplitude  observed  at  that  receiver.  Amplitude  is  nearly 
inversely  proportional  to  the  square  root  of  beam  density,  consistent  with  the 
calculation  of  the  geometric  spreading  by  ray  tube  area. 

All  combinations  of  sources  and  receivers  produce  ray  densities  that  are 
similar  in  form  to  Figure  7.  t.e..  uniform  over  broad  regions  surrounding  each 
receiver,  with  no  evidence  of  multipathing.  This  result  is  identical  to  that 
obtained  by  Comer  and  Aki  (1982).  They  found  that  multipaths  are  not 
generated  even  when  the  intensity  of  anomalies  in  the  Zandt  model  are 
doubled.  Multipathing.  however,  depends  on  the  scale  length  of  anomalies  as 
well  as  intensity.  This  is  emphasized  by  the  results  obtained  with  the  second 
r  del  investigated 


5.2  A  RANDOM  MODEL 


The  second  model  was  one  generated  by  perturbing  a  1-D  model  at  10  to  20 
km.  grid  points  in  horizontal  planes  with  a  rms  velocity  fluctuation  of  0.8°S 
(McLaughlin  and  Anderson.  1985)  Lnlike  the  Zandt  model,  this  model 
introduced  caustic  and  multipaths  at  teleseisrruc  range  This  made  it  essential 
to  use  Gaussian  beams  rather  than  asymptotic  ray  theory  (ART)  to  synthesize 
seismograms  at  receivers  in  the  vicinity  of  caustics.  ART  evaluates  the 
sup<  rposition  mt  eg?  i,  of  eq  t  lb)  by  a  stationary  or  saddle  point  lpproximat  ;or; 


The  stationary  pha- 
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LARGE  AMPLITUDE  STATION 


SMALL  AMPLITUDE  STATION 


Figure  7.  The  end  points  of  rays  at  the  surface  of  the  earth  at  70°  at  the 
highest  and  lowest  amplitude  stations  (azimuths  101°  and  281°  respec¬ 
tively)  in  the  experiment  shown  in  Figure  5.  The  1000  x  1000  km  region  of 
end  points  shown  heare  approximately  bounds  the  region  in  which  decay  of 
Gaussian  beams  is  less  than  e  '  at  0  03  Hz. 


tracing  problem  between  source  and  receiver,  leading  to  amplitudes 

proportional  to  factor  ===-.  This  factor  approaches  infinity  near  the 

vdetQ* 

caustic  surfaces  defined  by  detQ*  =  0  .  The  superposition  integral  and  its 
integrand,  however,  remain  regular  at  caustics  for  a  generalized,  complex  li 
matrix  (terven^  et  al.,  1982;  Cerven^,  1985a, b). 

5  2  1  Effect  of  varying  source  Location  at  constant  azimuth 

Figure  8  shows  contours  of  velocity  at  the  bottom  of  the  random  model  and 
the  location  of  sources  in  a  variable  source  experiment.  The  velocity  contours 
have  been  left  unlabeled  and  are  shown  only  to  illustrate  the  dramatically 
different  scale  lengths  of  velocity  fluctuation  in  this  model  compared  to  the 
Zandt  model.  At  a  constant  azimuth,  the  amplitude  fluctuations  (Figure  9)  due 
to  variations  in  source  location  are  both  larger  and  occur  more  rapidly  than 
those  seen  in  the  Zandt  model  (Figure  6).  This  reflects  the  fact  that  the 
smallest  scale  length  of  velocity  fluctuation  (10  km.)  roughly  equals  the  spacing 
of  source  points.  Greater  frequency  dependence  of  the  amplitudes  is  also  seen. 
This  is  due,  in  part,  to  the  presence  of  caustics  in  the  vicinity  of  the  receivers 
for  some  of  the  source-receiver  paths.  Figure  10  is  a  plot  of  ray  end  points, 
illustrating  the  development  of  one  of  these  caustics  in  the  vicinity  of  the  70° 
station  for  a  source  at  location  510-.  A  triplicated  zone  of  end  points  can  be 
seen,  which  is  elongated  along  a  narrow  azimuthal  zone.  Rays  having  end  points 
within  this  zone  are  found  to  have  a  one  unit  advance  in  their  KMAK  index, 
indicating  that  these  rays  have  passed  through  a  caustic  once.  A  receiver 
located  within  this  zone  of  triplicated  end  points  is  likely  to  observe  some  phase 
distortion  in  its  waveform  because  some  of  the  beams  that  contribute  to  the 

superposition  integral  will  have  a  phase  shift  This  phase  distortion  is 


Figure  8.  P  velocity  contours  in  a  horizontal  plane  at  the  bottom  of  a  model 
constructed  by  adding  a  0.8%  random  perturbation  to  the  velocities  of  a  1-D 
model.  Perturbations  were  assigned  at  10  to  20  km.  spaced  grid  points  in 
the  horizontal  plane.  The  contour  interval  is  0.05  km,  sec,  which  is  the 
same  as  that  shown  in  F  igure  4  The  epicentral  location  of  sources  530-  to 
530+  at  9  6  km.  depth  are  projected  onto  this  plane 


short,  period,  and  long  period  synthetic  P  waves  at  TO 


difficult  to  observe  in  synthetics  calculated  for  a  profile  of  stations  r400-  to 
r400+  shown  in  Figure  11.  The  phase  distortion  appears  as  small  change  in  the 
rise  time  of  the  broadband  pulse  at  station  rO.  The  small  negative  first  break, 
best  visible  on  the  broadband  pulse  at  rO,  is  not  due  to  the  phase  shifted  beams, 
but  rather  due  to  the  unique  interference  effects  of  this  particular  beam 
pattern.  The  largest  effects  to  be  observed  on  waveforms  might  be  expected  on 
the  short  period  instrument.  The  highest  frequency  band  would  have  the 
highest  per  cent  contribution  from  beams  within  the  triplicated  region  at  rO 
because  beams  more  distant  from  rO  would  suffer  a  stronger  exponential  decay. 
No  substantial  modifications,  however,  are  seen  in  the  short  period  waveform  of 
station  rO.  Comparison  of  amplitudes  in  different  frequency  bands  in  Figure  11 
now  shows  substantial  frequency  dependent  effects.  Long  period  amplitudes 
vary  only  about  half  as  much  as  short  period  and  broad  band  amplitudes.  A 
much  broader  area  of  beams  contribute  to  the  long  period  response,  smoothing 
over  the  effect  of  the  caustic  region  near  station  rO. 

Although  the  random  model  was  constructed  to  simulate  the  amplitude 

variations  seen  across  arrays  of  sources  or  receivers  having  apertures  of  the 

order  of  100  km.  (McLaughlin  and  Anderson,  1985),  it  does  not  predict  the  large 

variations  in  short  period  in  amplitudes  seen  over  much  smaller  receiver 

spacings  such  as  the  7  km  aperture  subarrays  of  NORSAR  (e.g.,  Thomson,  1983). 

A  random  layer,  having  the  same  statistics  and  thickness  as  that  placed 

beneath  the  source,  can  also  be  placed  beneath  the  receiver  array,  but  the 

combined  model  would  still  be  unable  to  reproduce  the  40%  variations  in 

amplitude  commonly  seen  across  the  NORSAR  subarrays.  Clearly,  smaller  scale 

lengths  of  velocity  fluctuation,  as  well  as  multiple  scattering  are  necessary  t.o 

explain  the  fluctuations  observed  over  smaller  aperture  arrays  Modeling  of  the 

effects  of  smaller  scale  lengths  was  not  possible  in  this  study  because  the 

asymptotic  and  paraxial  approximations  of  t  tie  ray  and  beam  techniques  begin 
to  breakdown  at  scale  lengths  on  the  order  of  or  less  than  the  wavelength 
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are  11.  Synthetic  P  waves  for  the  profile  of  stations  in  the  vicinity  of  the 


5  22  The  features  that  generate  caustics  and  multipaths 

Exactly  what  feature  of  the  random  model  was  responsible  for  the  thin, 
elongated  caustic  intersection  shown  in  Figure  10  ?  Since  this  feature  is 
elongated  along  a  particular  azimuth,  the  lateral  location  of  the  structure  is 
constrained  to  be  along  the  ray  paths  that  leave  the  3-D  portion  of  the  model  at 
this  azimuth.  The  range  of  vertical  take-off  angles  along  this  azimuth  is  also 
bounded  by  the  apparent  edges  of  the  caustic  intersections  in  the  plot  of  ray 
end  points.  The  calculation  of  the  KMAH  index  can  also  be  used  to  identify  the 
particular  rays  that  are  tangent  to  the  caustic  surface  at  depth.  By  either  of 
these  methods,  the  rays  that  describe  the  caustic  can  be  identified  and  their 
trajectories  plotted  through  the  3-D  region  of  the  model.  When  this  is  done 
(Figure  12).  it  can  be  seen  that  the  structure  responsible  for  the  caustic  at 
teleseismic  distance  is  a  low  velocity  zone,  extended  in  the  vertical  direction. 
The  reason  why  such  structures  have  been  generated  in  this  random  model  is 
that  the  grid  spacing  at  which  velocities  were  assigned  was  much  larger  in  the 
vertical  direction  (30  km.)  than  in  the  horizontal  direction  (10  or  20  km.)  Thus 
there  occasionally  will  be  regions  of  the  model  where  negative  perturbations 
strongly  correlate  between  adjacent  vertical  grid  lines,  forming  a  vertically, 
elongated  zone  of  low  velocities.  Similarly,  elongated  zones  of  high  velocities 
will  be  formed.  The  surprising  observation  seen  with  this  particualr  model  is 
that  only  a  very  small  perturbation  of  velocity  (0.8%),  with  a  vertical  scale  of  30 
to  60  km.,  and  a  horizontal  scale  of  10  km.  can  generate  caustics  and  phase 
advances  at  teleseismic  distances.  These  caustics  intersect  the  surface  of  the 
Earth  at  teleseismic  range,  but  the  areal  extent  of  these  intersections  are  too 
small  to  be  visibly  identified  except  for  some  rather  subtle  effects  in  the 
waveforms  at  a  few  number  of  stations.  The  generation  of  these  caustics 
depends  on  the  strength  of  velocity  fluctuation  as  well  as  on  the  relation 


teleseismic  range.  Contours  of  velocities  are  shown  at  a  0.05  km/sec  inter- 


between  the  characteristic  vertical  and  horizontal  scale  lengths  of  the  3-D 
modeL  In  the  example  considered,  the  shortest  scale  length  in  the  vertical 
direction  exceeds  that  in  the  horizontal  direction,  a  situation  which  is  probably 
not  the  common  state  of  crust/lithospheric  structure  (McLaughlin,  personal 
communication).  Some  notable  exceptions  to  this  would  include  regions  having 
concentrations  of  intrusive  pipes  and  plumes.  The  results  of  the  single 
modeling  experiment  described  here  suggest  that  some  distributions  of 
heterogeneity  would  produce  unacceptably  large  effects  on  teleseismic 
waveforms.  One  such  effect,  suggested  in  the  example  shown  in  Figure  11,  is 
that  the  first  motions  of  the  P  waves  from  explosive,  isotropic  sources  might  be 
interpreted  to  be  dilatations  rather  than  compressions  in  certain  narrow 
regions  surrounding  a  caustic.  The  fact  that  few,  if  any,  dilatations  have  been 
documented  in  the  P  waves  recorded  from  underground  explosions  may  be 
evidence  that  the  particular  structure  that  generates  these  caustics  is  a  highly 
uncommon  feature  of  the  Earth's  crust  and  upper  mantle.  From  this  example, 
it  is  clear  that  forward  modeling  of  the  effects  of  very  general  distributions  of 
heterogeneities  will  be  useful  in  defining  the  "heterospectrum"  of  the 
lithosphere.  (The  heterospectrum  is  a  term  adopted  by  Wu  (1986)  to  embrace 
both  the  the  strength  of  velocity  fluctuation  and  its  three-dimensional  spatial 
spectrum). 

5.3  VALIDITY  OF  THE  RESULTS 

All  of  the  example  seismograms  were  calculated  by  superposition  of 
Gaussian  beams.  The  accuracy  of  this  technique  depends  both  on  the  validity 
of  using  the  first  term  in  an  asymptotic  solution  to  the  elastodynamic  wave 
equation  and  making  a  Taylor  expansion  of  the  complex  phase  about  the 
central  ray  (the  paraxial  approximation).  It  is  thus  appropriate  to  question 


whether  the  3-D  models  of  velocity  discussed  in  this  paper  have  exceeded  the 
domains  of  validity  of  these  approximations.  Validity  constraints  in  continuous 
media  have  been  formulated  by  Beydoun  and  Ben-Menahem  (1985)  and  terven^ 
(1985b).  White  et  al.  (1988)  have  considered  problems  encountered  with 
continuous  media  as  well  with  boundary  interactions.  The  simplest  constraints 
to  check  are  those  related  to  the  asymptotic  expansion,  which  assumes 
decoupling  of  P  and  S  waves  to  zeroth  order  and  neglect  reflections  and 
conversions  by  regions  of  strong  gradient.  These  constraints  require  that 

wavelength  be  much  less  than  quantities  such  as  ~~ — -  and  — £- — -  (Kravtsov 
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and  Orlov,  1980;  Beydoun  and  Ben-Menahem.  1985),  where  v  and  p  are  velocity 
and  density  respectively,  and  i|  |  denotes  vector  length.  Both  the  Zandt  model 
and  the  random  model  satisfy  this  constraint  throughout  the  frequency  band 
0.03  to  4  Hz. 


Constraints  on  the  validity  of  the  paraxial  approximation  and  the  beam 
superposition  have  been  cast  in  terms  of  the  real  and  imaginary  parts  of  the  11 
matrix  and  distance  from  the  central  ray  (Beydoun  and  Ben-Menahem.  1985; 
Cerven^,  1985b).  In  order  to  see  if  these  constraints  are  obeyed  everywhere 
along  a  beam,  the  M  matrix  must  known  everywhere  along  the  central  ray 
associated  with  that  beam  At  any  point  0  along  the  central  ray,  11  (0)  can  be 
found  by  using  the  II  matrix  to  back  propagate  the  complex  11(0,)  matrix 
selected  at  the  end  point  00  Although  these  constraints  were  not  calculated  in 
the  two  examples  described,  a  check  on  the  overall  validity  of  the  beam 
superposition  was  taken  from  the  results  of  reciprocal  experiments,  in  which 
the  positions  of  sources  and  receivers  were  reversed.  Examples  of  this  test  in 
2-D  media  are  given  by  Nowack  and  Aki  (1984)  and  Muller  (1984)  Reciprocal 
experiments  were  conducted  in  the  form  of  allowing  a  plane  wave  to  be 
vertically  incident  on  the  3-D  models  The  plane  wave  was  expanded  into 


Gaussian  beams  by  the  procedure  described  by  Cerven^  (1982),  and  the 
wavefield  was  calculated  at  receivers  on  the  surface  of  the  3-D  models.  Such  an 
experiment  was  conducted  on  the  random  model  (McLaughlin  and  Anderson, 
1985)  and  on  a  model  having  a  heterospectrum  similar  in  scale  length  and 
velocity  fluctuation  to  the  Zandt  model  (Nowack  and  Cormier,  1985).  Although 
the  geometry  was  not  precisely  reciprocal  and  waveforms  were  not  directly 
compared,  both  experiments  produced  intensities  and  co-herence  of  amplitude 
fluctuations  similar  to  those  produced  by  the  teleseismic  experiments. 


Conclusions 
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The  theory  and  examples  discussed  in  this  paper  have  shown  how  the 
propagator  matrix  of  the  dynamic  ray  tracing  equations,  II ,  can  be  exploited  to 
connect  3-D  to  1-D  portions  of  a  model.  Doth  plane  wave  and  point  source 
initial  conditions  are  required  to  specify  the  elements  of  the  propagator  matrix. 
Thus  both  plane  wave  and  point  source  solutions  are  generally  useful  for  all 
asymptotic  methods  of  body  wave  synthesis  and  not  just  for  Gaussian  beams. 

An  investigation  was  made  of  the  effect  of  3-D  structure  in  the  source 
region  on  the  focussing  and  defocussmg  of  teleseismic  P  waves.  These  effects 
were  observed  with  two  different  models  as  a  function  of  source  location  within 
the  3-D  model  and  of  azimuth  at  the  source  In  a  3-D  model,  block  inverted 
from  teleseismic  travel  times,  ray  theoretical  amplitudes  matched  those 
predicted  from  superposition  of  Gaussian  beams  In  this  model,  the 
characteristic  scale  lengths  of  the  most  intense  velocity  fluctuations  (4%)  were 
on  the  order  of  50-100  km.  in  the  source  region  This  model  produced  a  factor 
of  two  fluctuation  in  amplitude,  associated  with  fluctuations  in  travel  time  on 
the  order  of  several  tenths  of  a  second.  Amplitude  variations  due  to  variations 
in  source  location  were  small  over  location  variations  that  were  small  with 
respect  to  the  scale  length  of  velocity  fluctuation  All  amplitude  variations  were 
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nearly  independent  of  frequency  across  the  body  wave  band  of  0.03  to  4  Hz. 

The  frequency  independence  of  amplitude  variations  across  the  body  wave  band 
may  have  important  implications  for  removing  the  effects  of  azimuthal 
amplitude  variations  due  to  3-D  structure  beneath  nuclear  test  sites.  Broad 
scale  length,  deep  seated  structure  can  affect  short  period  as  well  as  broad 
band  and  coda  measures  of  radiated  seismic  energy.  Its  effects,  however,  may 
be  easily  correctible  if  a  3-D  model  of  the  source  region  is  known  from  block 
inversion  of  travel  time  residuals.  A  resolvable  block  size  of  about  20  km.  may 
be  all  that  is  necessary  to  formulate  corrections  based  on  azimuth  of 
teleseismic  station  and  source  location  within  a  test  site. 

Calculations  with  a  random  velocity  mooel  demonstrated  that  a  smaller 
intensity  of  velocity  fluctuation  (0.8%)  can  produce  even  larger  amplitude 
variations  if  the  smallest  scale  length  of  fluctuation  is  on  the  order  of  10  km 
The  random  model  was  constructed  such  that  the  smallest  scale  length  of 
velocity  fluctuation  was  shorter  in  the  horizontal  direction  than  in  the  vertical 
direction.  \l  teleseismic  range,  this  model  generated  isolated  caustics,  which 
were  elongated  along  the  azimuth  of  approach.  The  waveform  distortion 
associated  with  these  caustics  was  small.  The  fact  that  the  caustics  were 
generated  at  all  by  such  mild  3-D  perturbations  is  significant.  It  suggests  that 
this  type  of  synthetic  modeling  may  be  useful  in  limiting  some  of  the  attributes 
of  the  heterospectrum  of  the  Earth's  lithosphere. 
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ABSTRACT 


The  focusing  and  defocusing  of  teleseismic  P  waves  predicted  by  a  known 
(Taylor,  1983)  three-dimensional  structure  beneath  the  Nevada  Test  Site  is 
calculated  by  dynamic  ray  tracing  and  superposition  of  Gaussian  beams.  The 
20  to  100  km.  scale  lengths  of  this  model,  having  velocity  fluctuations  of  several 
percent,  account  for  a  factor  of  three  fluctuation  in  the  azimuthal  patttem  of  P 
amplitudes.  Since  similar  sized  scale  lengths  and  intensities  of  velocity 
fluctuation  are  commonly  resolved  in  three-dimensional  inversions  of  mantle 
structure  in  other  regions,  focusing  and  defocusing  of  teleseismic  P  waves  is 
likely  to  be  a  ubiquitous  feature  of  every  test  site.  Hence,  network  averages  of 
that  weight  azimuthal  windows  in  which  P  energy  is  either  focused  or 
defocused  mil  tend  to  either  over  or  under  estimate  the  yields  of  underground 
tests.  The  results  obtained  with  a  known  NTS  structure  suggest  that  structural 
inversions  having  similar  structural  resolution  may  be  capable  of  reducing  25 
percent  of  the  variance  in  amplitudes  due  to  the  focusing  and  defocusing  of 
upper  mantle  structure  near  the  source. 

Assuming  that  scattering  processes  are  concentrated  in  the  crust  and 
uppermost  mantle  beneath  the  source  and  receiver,  the  integrated  energy  in 
the  P  coda  should  be  more  stable  than  direct  P  because  it  tends  to  homogenize 
the  effect  of  focusing/defocusing  structure  near  the  source  for  scattering  into 
direct  P  near  the  source.  For  direct  P  scattering  near  the  receiver,  however, 
even  the  late  coda  is  not  capable  of  completely  removing  the  effects  of 
focusing/' defocusing  near  the  source. 


INTRODUCTION 


The  underground  nuclear  tests  conducted  in  the  Pahute  Mesa  region  of  the 
Nevada  Test  Site  (Figure  1)  exhibit  a  unique  and  reproducible  azimuthal 
pattern  in  the  amplitudes  of  short  period  P  waves  at  teleseismic  range  (Lay  et 
al.,  1984).  Because  this  pattern  differs  from  those  of  tests  conducted  in  other 
regions  of  the  test  site,  it  cannot  simply  be  explained  by  either  receiver  effects 
or  by  variations  in  intrinsic  attenuation  along  different  paths  in  the  mantle. 

Evidence  of  tectonic  release  in  the  waveforms  of  long  period  P  and  S  waves 
of  Large  (greater  than  500  kt.)  Pahute  events  (Wallace  et  al.,  1983,  1984), 
suggests  that  tectonic  release  may  be  the  mechanism  that  also  produces  the 
azimuthal  variation  of  the  amplitude  of  short  period  P  waves.  A  study  by  Lay  et 
al.  (1984)  initially  favored  this  explanation,  but  noted  the  difficulty  in 
formulating  a  physically  realizable  mechanism  of  tectonic  release  that  would  be 
capable  of  strongly  affecting  short  period  amplitudes.  The  required  tectonic 
event  would  have  to  occur  nearly  simultaneously  with  the  explosion. 

Theoretical  and  nodel  studies  of  the  relaxation  of  the  likely  pre-stress  in  this 
region  have  found  negligible  effects  on  the  amplitudes  of  short  period  P  waves 
(Archambeau,  1972;  Bache,  1975). 

More  recent  studies  by  Lynnes  and  Lay  (1984;  1987)  and  Lay  and  Welc 
(1986a)  hypothesize  that  the  azimuthal  pattern  of  amplitudes  from  Pahute 
tests  is  instead  caused  by  focusing  and  defocusmg  of  P  waves  by  three- 
dimensional  structure  in  the  mantle  beneath  the  test  site.  The  structure  likely 
to  be  responsible  for  the  strongest  features  of  the  observed  focusing  and 
defocusmg  is  a  high  velocity  structure  in  the  mantle  beneath  Pahute  Mesa. 

Such  structure  was  first  proposed  by  Spence  (1974)  from  clusters  of  fast  P 
travel  times  from  Pahute  tests  to  teleseismic  stations.  Spence’s  model,  which 
qualitatively  fits  the  served  travel  time  residuals,  consists  of  a  conical  shaped 
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zone  of  positive  velocity  perturbation  of  6-8  percent  lying  directly  beneath  the 
Silent  Canyon  Caldera  (center  marked  with  x  in  Figure  1).  The  high  velocity 
perturbation  extends  to  190  km.  depth  and  its  half  width  is  about  50  km 

The  development  of  techniques  of  inverting  for  three  dimensional 
structure  (Aki  et  al.,  1977)  led  to  attempts  to  better  resolve  the  structure 
beneath  NTS.  including  the  the  structure  suggested  by  Spence's  study  (Minster 
et  al.,  1981;  Montfort  and  Evans,  1982;  Taylor,  1983).  A  consistent  result  of 
these  studies  is  that  the  high  velocity  structure  beneath  Pahute  Mesa  migrates 
to  the  north  as  depth  increases  (Figure  2).  The  strong  positive  perturbation  of 
velocities  in  this  structure  is  suggested  to  be  an  enriched  mafic  residue 
resulting  from  intermittent  melting  and  differentiation.  The  partial  melt 
material  itself  is  presumed  to  have  errupted  through  the  volcanic  center  of  the 
Silent  Canyon  Caldera. 

The  structure  resolved  in  the  inversions  of  Minster  et  al.  (1981)  and 
Montfort  and  Evans  (1982)  are  either  or  too  broad  in  scale  length  or  too  weak  in 
percent  velocity  fluctuation  to  have  much  of  an  effect  on  the  focusing  and 
defocusing  of  teleseismic  P  waves.  In  the  case  of  the  Minster  et  al.  model,  a 
smoothing  algorithm  reduces  velocity  perturbations  to  less  than  2  percent 
fluctuation  In  the  case  of  the  Montfort  and  Evans  model,  the  intent  was  to 
resolve  structure  over  a  much  broader  region,  using  larger  block  sizes.  Hence, 
the  approach  followed  by  Lynnes  and  Lay  (1984)  was  to  find  an  ad-hoc 
structure  consistent  with  the  observed  pattern  of  amplitudes.  The  resulting 
structure  bore  some  resemblance  to  a  hLgh  velocity  anomaly  mapped  in  the 
mantle  northeast  of  Pahute  Mesa  in  a  travel  time  inversion  by  Taylor  et  al. 
(1983)  The  effects  of  the  structure  resolved  in  the  Taylor  model,  however,  were 
not  directly  investigated 
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1.  A  representative  epicenter  of  an  underground  nuclear  test  at  Pahute 
Mesa  is  marked  by  the  cross  (x).  Iso-P  velocity  contours  ut  85  km. 
depth  are  shown  for  a  high  velocity  anomaly  resolved  in  the  study  by 
Taylor  (1983) 


Figure  2.  Iso-P  velocity  contours  greater  than  7.9  km./sec  for  the  Taylor  (1983) 
model  along  the  AA’  cross  section  shown  in  Figure  1.  Three- 
dirrensionally  traced  rays  (dashed)  are  projected  onto  the  cross  sec¬ 
tion  if  their  end  points  lie  within  20  km  of  the  plane  of  the  cross  sec¬ 
tion. 


This  paper  will  report  on  the  results  of  forward  modeling  of  teleseismic  P 
waves  for  an  explosive  point  source  placed  in  the  Taylor  (1983)  model  at  a 
Pahute  hypocenter.  For  several  reasons,  direct  calculation  of  the  focusing  and 
defocusing  effects  of  the  Taylor  model  is  relevant  to  reducing  the  uncertainties 
in  the  yield  estimates  of  underground  tests.  First,  the  model  would  be 
representative  of  the  best  information  that  would  likely  be  obtainable  at  foreign 
test  sites,  short  of  deploying  dense  local  arrays  within  the  test  site. 

Comparable  data  needed  to  obtain  similar  resolution  of  structure  would  consist 
of  a  regional  array  located  primarily  outside  the  test  site  and  a  well  distributed 
set  of  tests  with  known  origin  times  and  hypocenters.  Second,  the  effects  of  a 
known  model  for  the  structure  beneath  Pahute  Mesa  can  aid  in  identifying  any 
biasing  effects  of  focusing  and  defocusing  on  the  construction  of  empirical  yield 
verus  curves  from  NTS  data.  For  example,  if  globally  averaged  m^'s  are 
concentrated  in  a  narrow  azimuthal  band  of  focused  or  defocused  energy,  then 
yields  from  Pahute  events  will  provide  points  to  the  empirical  curve  that  are 
either  biased  toward  high  or  low  yield. 

EFFECTS  ON  DIRECT  P 

Seismograms  of  teleseismic  P  waves  were  synthesized  for  sources  in 
Taylor's  model  using  the  techniques  described  in  Cormier  (1986).  The  model 
was  specified  on  a  three  dimensional  grid  by  percent  perturbations  from  a 
reference  model.  Only  those  perturbations  were  used  whose  diagonal  elements 
of  the  resolution  matrix  exceeded  0.6.  Continuous  first  and  second  spatial 
derivatives  of  velocity  were  calculated  from  splines  under  tension  (Cline,  1981). 
Ray  tracing  was  performed  by  in‘<’grating  the  kinematic  and  dynamic  equations 
numerically  in  the  three-dimensional  region  surrounding  the  source  and 
analytically  along  the  rest  of  the  ray  path  m  a  flattened,  radially  symmetric 
earth  model  Seismograms  were  calculated  both  by  simple  ray  theory  and  by 
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superposition  of  Gaussian  beams  (terven^.  1985ab). 

Gaussian  beam  modeling  Superposition  of  Gaussian  beams  was  first 
performed  in  order  to  check  whether  the  three-dimensional  model  was  capable 
of  producing  any  frequency  dependent  effects  at  teleseismic  range.  Gaussian 
beams  remain  regular  at  caustics  and  can  succeed  in  predicting  the  first  order 
effects  of  frequency  dependence  at  caustics.  The  beam  superposition  also 
averages  the  effects  of  structural  scale  lengths  that  are  small  with  respect  to 
wavelength.  (Here,  structural  scale  length  is  taken  to  mean  the  distance 
between  local  minima  and  maxima  of  velocity  perturbation.)  Beam  width 
parameters  were  chosen  as  described  in  Cormier  (1986),  which  amounted  to 
taking  a  plane  wave  superposition  at  the  receivers  with  very  wide  Gaussian 
windowing. 

Broadband.  WWSSN  short  period,  and  WWSSN  long  period  seismograms  were 
synthesized  for  variable  azimuths  at  a  fixed  receiver  distance  of  60°  for  a 
Pahute  Mesa  source  located  as  shown  in  Pigure  1  The  synthetic  seismograms 
(Figure  3)  show  little  evidence  of  frequency  dependence,  with  the  relative 
amplitudes  of  waveforms  in  the  three  pass  bands  being  nearly  identical.  Similar 
results  were  found  in  a  study  of  a  three-dimensional  mode!  having  similar 
resolution  and  scale  lengths  of  structure  (Cormier,  1986)  Amplitudes  vary 
azimuthally  by  a  factor  of  three,  with  the  smallest  amplitude  at  due  north. 
Figures  1  and  2  show  that  ray  paths  in  northerly  directions  will  traverse  a  high 
velocity  anomaly  located  in  the  mantle  north  of  the  epicenter.  The  associated 
travel  time  variation  of  several  tenths  of  a  second  is  too  small  to  be  easily 
visible  in  the  synthetics.  Note  that  in  the  center  column  for  SP  WWSSN  in 
Figure  3,  the  reference  line  is  close  to  the  center  of  the  trough  of  the  waveform 
of  the  station  having  the  30  unit  amplitude  but  is  close  to  the  first  peak  of  the 
waveform  of  the  station  having  the  100  unit  amplitude  Ttie  sense  of  the  travel 
time  variation  is  consistent  with  focusing  and  defocusing:  small  amplitudes 
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correlate  with  fast  times;  large  amplitudes  with  slow  tunes  This  correlation  is 
discussed  and  shown  in  finer  detail  in  following  section  discussing  the  results  of 
ray  theoretical  modeling. 


Ray  theoretical  modeling  Since  the  Gaussian  beam  modeling  exhibited 
little  evidence  of  frequency  dependence,  it  was  deemed  appropriate  to  model 
seismograms  by  ray  theory  Ray  theoretical  amplitudes  were  predicted  from 
geometric  spreading  calculated  from  the  determinant  of  the  Q  matrix  obtained 
from  dynamic  ray  tracing  (Cerven£,  1985b).  A  dense  system  of  rays  was  shot  at 
variable  vertical  take-off  angle  and  azimuthal  angle.  The  range  of  vertical 
take-off  angles  was  chosen  to  correspond  to  great  circle  distances  between  35° 
and  95°.  The  determinant  of  the  Q  matrix  was  evaluated  at  end  points  of  the 
rays  at  the  earth's  surface.  The  ray  theoretical  amplitude  is  then  equal  to 

Figure  4  piots  these  amplitudes  versus  azimuth  for  the  system  of  rays 

shot  All  amplitudes  were  normalized  for  the  spreading  factor  appropriate  for  a 
P  wave  at  60°  in  a  radially  symmetric  earth  The  azim  .'h  plotted  is  the  azimuth 
at  the  point  at  which  a  ray  leaves  the  three-dimcns. ,;r.  '  region  surrounding  the 
source. 

Note  the  pronounced  minimum  in  ampht  .d*  -  .  ga:’i  occurs  at  northerly 
azimuths,  with  more  steeply  dipping  rays  at  gn  :  distances  having  the 
smallest  arnpht  udes  A  region  of  larger  amplit  u  !•  -  ■  urs  around  -70°  azimuth. 
In  other  azimuthal  ranges,  amplitudes  scatter  depend. ng  on  the  vertical  take¬ 
off  angle 

The  large  scatter  around  -  )0°  azimuth  (330  clockwise  from  north)  is  due 
to  the  presence  of  a  cau-t:  The  Gaussian  beam  synthesis  did  not  test  this 
azimuthal  range  Calculations  with  other  three-dimensional  models  obtained 
from  block  inversion  of  travel  times  and  with  statistically  generated  models  find 
that  such  caustics  can  be  easily  generated  by  veloc  ity  perturbations  as  small  as 
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one  percent  (Cormier,  1986;  McLaughlin  and  Anderson,  1987).  Tests  with 
seismograms  synthesized  by  superposition  of  Gaussian  beams  find  only  a  very 
small  amount  of  phase  distortion  in  the  P  waves  of  receivers  located  near  the 
caustic.  This  is  because  the  class  of  caustic  created  by  these  models  encloses  a 
very  small  area  on  the  earth's  surface.  This  area  is  generally  so  small  in  the 
vicinity  of  the  receiver  that  any  phase  distortion  or  frequency  dependence  in 
the  waveforms  is  barely  visible  within  the  pass  bands  of  common  seismographs. 

Comparison  of  predicted  with  observed  Amplitudes  Figure  5  shows 
observed  amplitudes  of  Pahute  events  collected  by  Lay  et  al.  (1984).  Note  the 
amplitude  minimum  predicted  by  the  Taylor  model  at  northerly  azimuths  is 
qualitatively  matched.  There  is  some  match  of  the  predicted  maximum  near 
-70°  azimuth  as  well  as  of  the  form  of  the  scatter  between  50°  to  150°  azimuth. 
The  minimum  in  the  observed  amplitudes  is  displaced  slightly  lo  the  east  of  the 
predicted  amplitudes.  It  appears  that  the  high  velocity  anomaly  in  the  Taylor 
model  is  placed  at  about  the  right  distance  from  the  center  of  Pahute  Mesa,  but 
is  displaced  to  the  west  relative  to  the  true  anomaly.  This  may  be  due  to  ray 
bending  in  the  vicinity  of  the  source,  which  the  inversion  technique  ignores. 
Even  without  any  azimuthal  shift  m  the  predicted  amplitudes,  the  predicted 
amplitudes  can  be  used  to  correct  for  focusing  and  defocusing,  giving  about  a 
25  percent  reduction  in  variance  u.  corrected  log  amplitudes.  It  is  interesting 
to  note  that  Taylor’s  model  predicts  a  similar  variance  reduction  in  travel  time 
anomalies.  This  suggests  that  a  particular  model's  success  in  reducing  the 
variance  in  travel  times  can  be  taken  as  a  rough  measure  of  its  potential 
success  ir,  reducing  the  variance  m  log  amplitudes  due  to  focusing  and 
defocusing 

The  strong  concentration  of  observations  around  the  amplitude  minimum 
will  tend  to  reduce  the  network  averaged  of  Pahute  events  This  effect  can 
be  alleviat  ed  by  more  sophist  seated  averaging  of  the  tndiv:  Mai  iris's  reported  by 


<s> 


© 

O) 

I 


© 

oo 


aanindwv 


.'.N  -V 


:.y.  -v 


« v  *  -.  . 
i"1'  I*~i^i1h  *  *h~  > 


AZIMUTH  (deg 


Crfif 


stations  Ln  a  global  network.  For  example,  maximum  likelihood  processing  can 
take  into  account  variations  in  the  azimuthal  distribution  of  stations  (Ringdal, 
1976)  and  hence  can  be  designed  to  provide  azimuthally  dependent  weights 
based  on  the  focusing  and  defocusing  properties  of  the  source  region. 

The  dependence  of  focusing  and  defocusing  effects  on  vertical  as  well  as 
azimuthal  ray  angle  is  better  illustrated  in  a  plot  of  amplitudes  on  a  focal 
sphere.  This  is  shown  in  Figure  6.  In  this  plot,  the  agreement  between  observed 
and  predicted  amplitudes  is  more  compelling,  with  some  second  order  features 
that  depend  on  both  take-off  angles  highlighted.  Note  the  agreement  of 
negative  anomalies  near  azimuths  75°.  120°,  150°,  and  190°  as  well  as  the 
relative  size  of  negative  and  positive  anomalies  in  the  northeast  quadrant. 

Correlation  of  predicted  travel  times  and  amplitudes  Figure  7  shows  the 
azimuthal  and  take-off  angle  dependence  of  travel  time  residual  predicted  by 
Taylor’s  model.  The  same  results  are  shown  in  the  form  of  foe  !  sphere  plots  in 
Taylor's  paper.  The  residual  is  taken  from  the  predicted  tiny,  of  PREM 
(Dziewonski  and  Anderson,  1981).  Figure  9  plots  this  travel  time  residual 
against  amplitude.  The  general  sense  of  the  correlation  expected  for  focusing 
and  defocusing  is  apparent,  with  slow  times  correlating,  .'.it!,  large  amplitudes. 
The  scatter  from  this  trend,  however,  is  large.  Thus,  w .  k  or  non-existent 

correlation  of  travel  time  anomalies  with  amplitude  a: . ulies  of  data  should 

not  necessarily  be  taken  as  proof  that  focusing  and  defocusing  is  not  affecting 
amplitudes.  Similar  weak  correlations  between  amplitudes  and  travel  times 
have  been  found  in  other  studies  (e  g.,  MeCreery,  1986). 

An  intuitive  argument  why  log  amplitudes  do  not  precisely  correlate  with 
travel  times  can  be  made  by  considering  travel  time  curves  and  ray  density 
plots  of  test  structures.  It  can  be  seen  that  the  common  situation  is  that  rays 
will  tend  to  cluster  near  causti'-  surfaces  and  regions  of  high  velocity  gradient. 
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Correlation  of  calculated  amplitudes  and  travel  times  for  teleseismic 
1  waves  from  a  Pahute  Mesa  source  in  the  Taylor  model. 


The  rays  having  the  slowest  travel  times,  however,  will  tend  to  be  slightly 
displaced  from  the  regions  of  densest  rays.  An  example  of  this  can  be  seen  in 
the  fault  ione  calculations  of  Cormier  and  Spudich  (1933).  In  that  example,  the 
highest  amplitudes  are  associated  with  caustics  near  a  region  of  high  lateral 
gradient  at  the  edges  of  the  fault  zone  but  the  slowest  travel  times  are  in  the 
center  of  the  fault  zone.  Thus,  log  amplitudes  will  tend  to  be  correlated  only  in 
a  general  regional  sense  with  travel  times  A  travel  time-log  amplitude  plot  will 
typically  have  large  scatter. 

Effects  of  model  parameterization  Velocity  models  parameterized  by 
blocks,  while  convenient  for  travel  time  calculations,  are  not  the  most  suitable 
for  amplitude  calculations.  The  sharp  velocity  contrasts  at  block  boundaries 
generate  clusters  of  low  and  high  amplitudes  that  rapidly  vary  over  the  focal 
sphere.  Moreover,  sharp  boundaries  at  the  edges  of  bl”"ks  ::  -  nig  linear 
dimensions  on  the  order  of  10  to  20  km.  may  not  be  repress  .  ‘  dive  of  the  true 
variation  of  earth  structure  at  these  scale  lengths  m  the  up;  r  mantle  For 
these  reasons,  the  Taylor  model  was  instead  par  arm  ter. zed  L>  -plim  s  under 
tension  (Cline,  1991).  P  velocities  were  specific  1  knots  placed  at  the  center 
of  the  original  bloi  k'  When  the  tension  paraclete :  J  -  0,  velocities  are 
interpolated  by  mb.  polynomials  When  o  -  32  .  th<  interpolator]  is  essentially 
linear 

Because  ray  paths  are  strongly  perturbed  by  changes  m  velocity  gradient 
and  the  dynamic  equations  governing  geometric  spreading  depend  on  the 
second  spatial  derivative  of  velocity,  one  may  expect  th  it  the  tension 
parameter  would  uTe'-t  the  t*de$etsmie  amplitude  pattern.  A  :>  <1  of  ttus  is 
shown  :n  Figure  q.  whicr.  •  ornpar**s  predicted  arrplitudes  for  a  =  10  0  and 
(7  =  00  Ihe  broad  regions  cf  negative  and  positive  perturbations  remain 
substantially  the  same  m  the  t  wo  eases  Negative  anomalies  are  observed  in 
botti  cases  to  the  north  and  northeast,  with  more  negative  anomalies  observed 
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for  more  distant  stations  at  steeper  vertical  take-ofT  angles.  The  most 
pronounced  differences  are  seen  to  be  associated  with  caustics,  marked  by 
narrow  regions  of  tightly  dusted,  positive  anomalies.  For  example,  note  the 
differences  around  the  azimuths  bounded  by  300°  and  340°.  !n  the  region  of 
the  model  producing  the  anomaly  pattern  in  this  zone,  zero  tension  introduces 
a  high  velocity  region  between  two  different  low  velocity  regions.  The  amplitude 
anomalies  between  the  two  caustics  become  weakly  negative  instead  of  weakly 
positive  In  essence,  low  values  of  tension  tend  to  decrease  the  fundamental 
scale  length  of  the  model  by  allowing  strong  fluctuations  of  velocity  between 
knot  points.  Clearly,  a  higher  tension,  corresponding  to  linear  interpolation  is 
desirable  because  it  will  tend  to  be  more  faithful  to  the  original  scale  lengths 
resolved  by  the  travel  ♦  .me  inversion 

Simple  ray  theory  rather  than  superposition  of  Go  issian  beams  was  chosen 
for  the  calculations  m  F.gures  6  and  9  because  it  enabled  a  r  ■  ■  rapid 
calculation  of  the  amplitude  anomalies  at  a  high  der.sdy  of  «•  •  ins  and 
because  there  was  no  i  videm-e  of  strong  frequency  d*  pende*  •  n  a  trial 
synthesis  by  bean:  -uperpo^ition  If  the  amplitude  ammalies  shown  in  Figure  9 
had  bee:,  tl  : '  .  ...  „■  mi  superposition  rat!  •  •  _>  i ‘  heory .  an  even 

closer  comp!',  ,  .  ;..r.e  beer:  obtained  i  ms  -  1  .  ;?<  superposition 

tends  t  ■  iv.  r  ;g*  ?  tit  s* ;  .  tore  over  a  Fresnel  /•>".<  w‘  :  .  fo:  t :  <•  short  period 
band,  is  on  tie.  order  of  U.e  distance  between  knot'  m  the  beam 

superposition,  the  tight  clusters  of  posit. vr  arnphtud'  nioniah'  -  seen  near 
azimuths  300°  arid  340°  wouid  be  reduced  iri  r  1  ~pr<  ad  o.,t  over  a  slightly 

i 

larger  reg.or.  As  d,se;. •-.<.••  <•  lr'.ier,  these  im.,0  .--s  ar*-  *  io  sp,it..:l.y 

com-ent  rut  ml  '  >  prod  ;•••  my  phase  distort.  »r.  m  'to  ;  i<s  band  of  short  period 

.  instruments 
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Coda  stability  and  waveform  complexity  For  variations  in  source  and 
receiver  site  on  the  order  of  10  to  100  km,  magnitude  measures  based  on 
integrated  P  coda  vary  less  than  those  based  on  the  amplitude  of  the  first 
several  cycles  of  the  P  wave  (Ringdal,  1983;  Baumgardt,  1985;  Gupta  et  al., 

1986).  This  observation  is  consistent  with  the  effects  of  focusing  and 
defocusing  by  large  scale  structures  in  the  mantle  beneath  the  source  region 
together  with  the  assumption  that  the  coda  is  generated  by  scattering  by 
heterogenities  concentrated  in  the  crust  and  uppermost  mantle.  In  the  source 
region,  the  effects  of  focusing  and  defocusing  are  minimized  for  scatterers 
located  at  a  distance  from  the  source  that  is  greater  then  the  characteristic 
wavelength  of  heterogeneity  in  the  mantle  beneath  the  source  (Figure  10).  For 
a  hypothesis  of  the  c  oia  being  dominated  by  single  scattering,  the  later  part  of 
the  coda  will  become  progressively  better  m  minimizing  the  effects  of  focusing 
and  defocusing  This  is  because  the  later  coda,  being  general  -  d  by  scattering 
structure  further  fr  t*v  epicenter,  is  more  likely  to  travel  to  the  receiver  by 
a  P  wave  path  th  r  •*•••  -  sample  the  same  man  lie  heterogeneity.  These 
different  r i u » : . t ! •  he!  t  !,.•<  will  produce  differ*  tit  d*  grees  of  of  focusing  and 

d«'fo"  is.ng  1 1  ,i  •  ’  .  '  heterogeneity  samp!*  !  by  th'  P  w.c.c  path 

r*'  -  p',:.  -  .hie  ’  ■’  ‘  >  •  r  ■  ■.  a!  In  thro  tit:,*  :  -  .  g*.  on.,  try.  this  effect  is 

magti.fi'  •!  b>  •  ••  v.ittcruijj  th  it  ot  vwthm  a  concentric  ring 

surrounding  the  sour: 

A  simplified  cal'1  ui.it ; on  of  this  effect  for  the  Pahute  structure  is  shown  in 
Figures  1  1  and  12  The  s<  at  ter.ng  is  assumed  to  be  doriun.it  id  by  scattering  of 
energ,  prop. ig.it  .rig  ar  •  -.  ir./mit-illv  . n  the  ru-l  us  c.  Ig  or  Kg  waves  into  P 
wave.  pr  j  ..g  it  -,g  -  •  v  iw.iy  fr  im  the  source  region  A  ring  of 

sc  it  t  er,.rs  i--,,rr.>  !  it  I'.idu.-’  if  D  km  from  a  Puhutr  sour'  e  1’hc  radius  of 
the  ring  is  <  t.oser;  to  - .  r  r : .  i : .  i  *  *  •  ‘tv  <  ffe<  t  -  on  the  o  la  arriving  19  sec  after  the 
direct  P  wave  f  ir  !  g  »r  kg  ’  ■>  P  at  termg  The  group  velocity  of  t  he  S.  Ig,  or 
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Figure  10  Schematic  ray  paths  of  direct  P  and  P  r-oda  energy  for  a  source 

located  near  a  velocity  heterogeneity  in  the  upper  mantle.  P,  S.  Lg.  or 
Hg  energy  propagates  in  the  direction  of  the  dashed  line  away  from 
t  e  source  and  is  scattered  into  a  P  wave  by  a  heterogeneity  in  the 
crust  or  upper  mantle.  J 
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Figure  11.  An  approximate  calculation  of  focusing  and  defocusing  effects  is 
made  for  coda  generated  near  the  source  by  considering  a  ring  of 
scatterers  located  33  km  laterally  from  the  source. 
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Rg  wave  is  assumed  to  be  3.3  km/sec.  Wide  angle  scattering  is  assumed  in  the 
vertical  angle  but  only  narrow  angle  scattering  is  assumed  ir;  the  azimuthal 
angle  Scattering  near  the  receiver,  scattering  m  the  deep  mantle,  and  multiple 
scattering,  are,  for  the  time  being,  neglected 

There  is  still  a  suggestion  of  the  northerly  minimum  in  coda  amplitude  in 
Figure  12,  but  the  minimum  is  not  as  intense  or  as  wide  in  azimuthal  range  as 
for  the  direct  P  wave  The  nearly  constant  levels  of  amplitude  in  westerly 
azimuths  for  more  distant  stations  are  due  to  the  fact  that  the  scattered  P  rays 
traversed  a  portion  of  the  model  in  which  the  reference  one- dimensional  model 
of  structure  was  unperturbed 

These  effects  are  consistent  with  those  observed  by  Lay  and  Welc  (l988ab) 
in  Pahutc  waveforms  in  the  3  to  10  sec  time  window  following  the  direct  P. 

They  found  the  coda  to  be  less  affected  by  defyus.ng  i r :  northerly  azimuths. 

The  effects  decreased  with  later  time  in  the  coda  P  wavefor:  s  were  generally 
more  complex  in  northerly  azimuths  Phis  .-gr,  ....  v  •  h  on**  of  the  scenarios 
proposed  for  the  origin  of  waveform  compb’x  t  by  *'  uglas  ct  al  (1973),  in 
which  complex  waveforms  ar**  created  when  tic  !.■  '  P  *  we  is  defocused  by 

rnant  i  e  structure  relative  to  t  he  P  cola 

Other  examples  >f  forward  m.uliling  f  .  :  0  ir:d  defoeusing  (Cormier, 

1986,  Mr ’aughlin,  1987)  find  that  variations  r.  -  >  location  .is  small  as  30  to 
f>0  km  are  sufficient  to  strongly  affect  t  fie  amj  dud-  s  of  direct  P  waves  The 
controlling  factor  is  the  characterist  ic  scab*  n  ngths  of  ]  percent  and  greater 
velocity  heterogeneities  in  ttic  upper  mantle  Tties*  observations  hear  upon  a 
hypot  tl»*s;s  o*  Douglas  «•»  ai  |-!1  96  1  .!  !'•>  .v  "o.n  ’  f  >r  ompb-x  and  simple 

waveforms  ,tjsrTV'*d  ovr  m  arly  id*  rite  al  rav  j  attis  Douglas  ct  at  suggested 
that  the  coda  was  genera!  rd  by  weak  scattering  ail  along  t  tie  ray  path  rattier 
than  by  scattering  cote  •Titrated  ,n  the  rust  and  uppermost  mantle  The 


-  A  fo'  di  sphere  plot  *>f  the-  amplitude*  of  rodrt  arriving  10  sec'  after 
dir^r't  P  using  th**  tppr  •'Ximat  »•  *x  tl»  o  [at  ion  illustrated  n:  Fig  urp  ’  1 
The  ring  of  si  VIitit*  <*as  ass  coed  *n  be  centered  around  the  epi¬ 
center  <howr.  .r.  V.H-..T  •  '  "’yirit,  .ire  defined  , r:  K^jurr  f>  This  r.ilru 
lation  only  simulates  effec  ts  on  <  oda  generated  near  t  he  source  The 
coda  will  also  have  a  ■  rnpouent  generated  near  the  receiver 


results  of  the  focusing  and  defocusing  experiments  reported  here  and  in  other 
studies,  however,  suggest  that  it  may  not  be  nec  essary  to  invoke  significant 
scattering  in  the  deep  mantle  to  account  for  simple  and  complex  waveforms  for 
nearby  ray  paths,  if  "nearby '  is  taken  to  be  distances  as  large  as  30  to  50  km 

Wide  versus  narrou.  angle  scattering  Wide  angle  scattering  in  the 
azimuthal  angle  can  be  more  effective  in  reducing  the  effects  of  focusing  and 
defocusing  compared  to  narrow  angle  scattering  Wide  angle  scattering  allows 
greater  opportunities  for  scattered  F  waves  to  sample  regions  outside  of  a 
particular  focusing,  defocusing  structure  in  the  upper  mantle  From  amplitude 
and  phase  correlations  at  NORSAR.  Wu  and  Flattc'  (1986)  suggest  that  wide 
angle  scattering  is  particularly  important  in  the  upper  1  to  krn  of  the  crust 
They  postulate  that  wide  angle  scattering  at  sha!!-*  depths  an  account  for 
imperceptible  fluctuations  m  travel  lime  being  as-o  lated  wi(v  ,  irge  amplitude 
fluctuations  at  small  (7  krn.)  aperture  sub  arra\  >•  (.•  g  Thorn-  1903) 
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Tins  ;s  because  the  small  region  in  which  rays  arc  focused  or  defocused  near 
the  source  is  spread  out  over  a  very  broad  region  after  propagation  to 
tclcseismic  range  Scattering  -an  occur  over  a  very  broad  region  near  the 
re  eiver  for  rays  having  a  very  small  distribution  of  vertical  and  azimuthal 
take-ofT  angles  of  direct  P  near  the  source.  Assuming  that  scattering  is 
uniformly  concentrated  m  the  crust  and  upper  mantle  near  the  earth's  surface, 
the  P  coda  for  an  explosion  or  earthquake  near  the  earth’s  surface  will  be 
composed  of  roughly  equal  amounts  of  P  energy  scattered  near  the  source  and 
near  the  receiver  The  P  coda  of  a  deep  focus  earthquake  will  primarily  be 
composed  of  P  energy  scattered  near  the  receiver.  These  predictions  have 
been  verified  by  comparison  of  the  phase  velocities  of  late  P  coda  energy  from 
deep  and  shallow  earthquake  sources  (Dainty.  1996).  With  this  model  of 
scattering,  it  is  clear  wh>  eo.la  magnitudes  of  shallow  ■  nts  will  tend  to  be  only 
about,  ffty  percent  successful  in  homogenizing  the;  effete.  of  focusing  and 
defocusmg  structures  near  the-  source  and  receiver.  At  each  receiver,  a  coda 
magrutud*  will  be  successful  ,n  removing  the  effects  of  'erasing  'defocusing 
struct  ares  beneath  the  receive  r,  with  th»‘  perform, un  •  r:  ogressively  improving 
with  .n  rr  ismg  t;me  into  t  he  oda  Even  tha  very  ’  >!•  j  .or,  of  the  coda, 
h-iwev.  r  *,•_>  ;.!)t  totally  -u  >.  •  m  g  tic  t  "  is  of 

fo.  ..sc  g  tef  irus.ng  "t  ru-  tore  t  ■  n*  ,»th  the  so  c  >•  This  is  because  about  fifty 
P*  r»-e-  ■  '  ate  i  oda  will  v  it.s.s!  of  d.t  'I  :  rnit  ha<  been  focused  or 
def  i  .  !  near  t h "  source  »•  1  '  t  ■  c .  -■  '  ■  •  ■  •  ;  ,,r  t  ui ver 
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more  stable  than  direct  P,  but  not  necessarily  for  every  source  and  receiver 
combination  (Bullitt  and  Cormier,  1984).  The  increased  stability  will  depend  on 
the  relative  importance  of  scattering  and  focusing/defocusing  near  the  source 
versus  that  occurring  near  the  receiver.  The  worst  performance  of  coda 
magnitudes  will  occur  in  a  situation  in  which  the  strongest  focusing/defocusing 
is  in  the  upper  mantle  close  to  the  source  and  in  which  scattering  in  the  crust 
near  the  source  is  relatively  weak. 

CONCLUSIONS 

A  known  three-dimensional  structure  beneath  Pahute  Mesa,  Nevada  Test 
Site  (Taylor,  1983)  can  account  for  many  of  the  features  in  the  azimuthal 
amplitude  pattern  of  teleseismic  P  waves  from  Pahute  underground  tests.  This 
model  can  be  used  to  correct  for  focusing  and  defoeusing  effects  of  the 
structure  beneath  Pahute  Mesa  accounting  for  factors  of  three  in  amplitude 
fluctuation  and  for  0.6  sec.  in  travel  tune  fluctuation.  "rhe  reduction  in  variance 
of  teleseismic  magnitude  or  log  amplitude  using  these  corrections  is  about  25 
percent,  similar  to  the  reduction  of  variance  in  teleseismic  travel  times.  These 


results  are  useful  in  predicting  the  structural  resu’.-tic::  needed  for  models  of 
other  test  sites  to  be  useful  in  making  corrections  for  r  using  and  defocusing. 
The  NTS  results  suggest  that  meaningful  co: -"ct.io;. •  •  o  be  made  if  the  model 
resolves  10  to  20  km  scale  lengths  down  to  ’00  km  will,  perturbation^  of 
velocity  exceeding  4  percent.  Velocity  inversions  that  primarily  resolve  scale 
lengths  larger  than  these  or  that  smooth  over  anomalies  larger  than  2  percent, 


tests  spaced  less  than  10  km.  apart.  It  is  also  necessary  to  obtain  an  average 
crustal  structure  within  the  test  site  from  seismograms  recorded  at  local  and 
regional  range. 


The  focusing  and  defocusing  effects  of  20  to  50  km.  scale  length  structure 
having  perturbation  in  P  velocity  of  several  percent  is  nearly  independent  in 
the  frequency  band  of  teleseismic  body  waves.  This  conclusion  is  even  stronger 
in  the  0.2  to  10  Hz.  band  in  which  measurements  are  made  on  the  teleseismic  P 
waves  of  underground  nuclear  tests.  Frequency  dependent  effects  in  the  coda 
of  teleseismic  P  waves  are  probably  due  to  either  the  effects  of  heterogeneities 
having  scale  lengths  smaller  than  20  km.  and/or  to  frequency  dependent 
effects  in  the  scattering  processes  occurring  near  the  source  and  receiver. 

Magnitude  measurements  based  on  the  integrated  energy  in  the  coda  of 
teleseismic  P  waves  are  likely  to  be  more  stable  because  they  can  remove  some 
of  the  focusing/defocusing  effects  of  three-dimensional  mantle  structure  near 
the  source  The  deeper  in  the  coda,  the  measurement  is  made,  the  less  affected 
it  will  be  by  mantle  structure  near  the  source.  The  optimal  time  in  the  coda  for 
this  measurement  should  be  as  long  as  possible  after  the  direct  P  wave  given 
the  signal  to  noise  ratio.  The  minimum  time  to  achieve  good  stability  can  be 
estimated  by  dividing  the  length  of  characteristic  scale  lengths  of  mantle 
structure  near  the  source  by  the  velocity  of  the  presumed  scattered  wave  near 
the  source.  For  example,  assuming  a  3.3  km./sec.  S  wave  is  scattered  into  a  P 
wave  that  is  propagated  to  teleseismic  range,  one  would  est.mate  that  after  30 
sec.  into  the  coda,  the  focusing/defocusing  effects  of  100  km.  scale  and  smaller 
length  structure  beneath  the  source  would  be  minimiz'  d  Coda  magnitudes, 
however,  are  only  partially  successful  in  removing  the  focusing  defocusing 
effects  of  structure  beneath  the  source  rhey  cannot  remove  '  hrse  effects  from 
the  friction  of  the  coda  that  is  due  to  scattering  of  direct  P  m  ar  the  receiver. 


These  conclusions  are  consistent  with  tests  of  the  relative  performance  of 
coda  versus  classical  magnitudes.  The  predicted  behavior  of  coda  amplitude 
critically  depends  on  assumptions  about  the  distribution  of  mantle 
heterogeneity  with  depth.  Smaller  scale  heterogeneities  with  greater  percent 
velocity  fluctuations  are  assumed  to  be  concentrated  closer  to  the  surface. 
Scale  lengths  on  the  order  of  several  kilometers  to  10  kilometers  are  assumed 
in  the  crust  and  scale  lengths  on  ther  order  of  10  to  100  kilometers  are 
assumed  in  the  upper  mantle.  Fluctuations  in  the  mid  and  lower  mantle  down 
to  the  D"  layer  near  the  core  are  assumed  to  be  smaller  than  1  percent. 
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t>  Three-dimensional  slab  effects  on  S  waves 


5.2  A  FEASIBILITY  STUDY  USING  VERY  WIDE  GAUSSIAN  BEAMS 

It  is  quite  resonable  to  assume  that  focusing  and  defocusing  by  slabs  will  be  a 
frequency  dependent  phenomena.  The  slab  is  a  narrow  structure  and  one  would  expect 
that  high  frequency,  short  wavelength  radiation  would  feel  the  presence  of  the  slab, 
whereas,  low  frequency,  long  wavelength  radiation  would  average  out  or  not  see  the 
effects  of  a  slab  wider  than  the  wavelength.  To  calculate  these  effects,  it  is  necessary  to 
use  a  technique  which  does  not  use  the  ray-theoretical  stationary  phase  approximation. 
Instead,  a  technique  is  needed  that  incorporates  waves  arriving  at  many  different 
directions  and  frequencies.  Gaussian  beam  superposition  and  WKBJ  ' Maslov 
seismograms  belong  to  a  class  of  asymptotically  approximate  techniques  that  can 
include  these  effects  Comprehensive  reviews  of  these  respective  techniues  are  given  in 
f'erven^’  (1985a.b)  and  Chapman  and  Drummond  (1982)  and  Thomson  and  Chapman 
(  19HR)  Comput  ational  details  associated  with  the  examples  shown  here  are  described  m 
s  ibsi'  -t  ion  5  2 

Figure  5  1  s flows  the  results  of  synthesizing  an  S  wave  at  50  great  curie  degrees 
awav  from  i  "'()()  kiii  deep  earthquake  in  i  f>()  d*  gree  dipping  slab,  which  penetrates  to 
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Synth^'ic  oF  at  50  great  circle  degrees  and  variable  anmuths 

from  a  ?00  <m.  deep  earthquake  located  m  a  60  degree  dipping  slab 
that  penetrates  to  1200  km.  flab  parameters  are  taken  from  Creager 
and  Jordan  ( 19H6)  for  the  Kunl-Kamchatka  slab.  Computational  detail 
of  the  synthesis  are  described  in  Appendix  A  of  this  proposal  and  in 
Cormier  ( 1965b)  Fhe  numbers  at  the  left  of  each  pulse  measure  peak 
amplitude  M  ■  r 
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Several  slab  effects  are  illustrated  in  Figure  b  1  First.  stations  at  a/miutlis  on  the 
dipping  side  of  the  slab  (<?7l)°  and  arrive  earln-r  than  stations  on  tt:e  side  away 

from  the  dipping  direc  tion  (4‘o°  and  90°)  Second,  station  at  a/unut fis  on  the  dipping 
side  have  been  defoeussed  by  the  high  velocity  slat;  and  art  -mailer  arnpld  ude  compared 
to  stations  at  azimuths  on  '.tie  side  away  from  the  dipping  d.reetion  (the  numbers  77,  H  i 
et  e  .  indicate  peak  relative  amplitude)  Third,  the  stations  on  t  tie  clipping  side  exhibit  a 
tan  >f  energy  that  persists  for  up  to  10  or  more  seconds  after  the  peak  1‘tns  tail  rr.av  tie 
termed  a  slab  diffract  ed  wave  At  the  station  at  ')Ti"  a/.muth.  the  tail  <  xibits  a  a 
so  ondary  peak  quite  see  i  r  to  t  tie  examples  of  slab  n'Tr  a.  1  .or.  shown  .  r .  t  fie  two 
dimensional,  finite  difference  calculations  of  Vidale  (19B(>) 
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:  ;  .  ■■  -..mm  :  over  vertical  ami  azimuthal  take  off  antles.  r  is  simply 

V  i'.  <M\:f  leiay  time  modified  by  an  imaginary  perturbation  that  depends  on  the 

*  'ii.-i::- s:  in  beams  and  the  distance  to  the  station  measured  in  ray 

c  ordinal es  An  additional  damping  factor  has  been  included  (Madariaga, 
ei;  M.cda'miga  and  Papadirnt  riou,  199b)  that  depends  on  the  sampling  rate  A<  .  The 
••  it r  M.  Q,  -lr  arc  calculated  by  dynamic  ray  tracing,  which  involves  integration  of 
i  sc-  if  ,me.ir  equations  along  a  ray  path.  These  equations  are  in  addition  to  the 
m  ncriiat ,  e  equations  for  ray  trajectory  and  S  polarization  (p.g..  Clcrven^  and  I  iron,  19Hi> 
billies  1994;  fmrveny,  1985a,b)  The  real  part  of  the  matrix  H  was  chosen  to  corrca  f  v 
the  effects  of  gradients  in  the  receiver  region  and  makes  the  Gaussian  beam 
superposition  essentially  equivalent  to  WKBJ /Maslov  superposition  of  plane  w-iv  - 
rec  eiver  witti  a  small  amount  of  Gaussian  windowing.  The  quantities  /■,  c 
components  of  vector  slowness  evaluated  at  the  receiver  The  m-.-*  •  • 
problems  in  computation  are  to  make  sure  t  hat  a  dense  >  ,g>  ■  • 
that  the  complex  delay  time  r  is  sampled  at  !*-<*>:  it  •  *.•  •  «•  • 

The  sampling  problems  are  quite  analguus  t  .  -•  :■ 

(1979)  for  V\KB!  seismograms  m  a  v- ■»-*  , 

be  interpolated  in  amplitude  A  e.  :  - 
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angles  at  which  the  discrete  superposition  is  terminated.  In  a  3-D  synthesis  there  are  8 
types  of  truncation  phases,  corresponding  to  4  sides  and  4  corners  of  the  region  covered 
by  ray  end  points.  Narrow  beam  widths  can  remove  these  truncation  phases,  but  can 
also  sometimes  remove  physical  radiation.  A  known  example  is  that  of  head  waves  rays, 
in  which  a  broad  spectrum  of  angles  is  required  to  properly  describe  the  head  wave. 
Here,  as  shown  in  Figure  5.4,  a  so  called  "optimal”  beam  width  (Klimes*v,  1986)  removes 
physical  radiation  corresponding  to  slab  diffraction.  The  consensus  among  investigators 
looking  at  a  variety  of  2-  and  3-D  problems  is  that  in  order  to  include  phenomena 
requiring  a  broad  angular  spectrum,  one  must  use  as  wide  a  beam  width  as  feasible  and 
make  the  technique  approach  the  WKBJ/Maslov  limit  of  plane  wave  superposition 
(Felsen,  1984;  Madariaga,  1984;  Wu,  1985;  Thomson  and  Chapman,  1986;  Lu  et  al.,  1987). 
Truncation  phases  can  to  some  extent  be  reduced  by  a  small  amount  of  Gaussian  beams, 
or  as  suggested  by  Thomson  and  Chapman  (1986)  be  minimized  by  asymptotically 
estimating  the  integrals  in  the  total  angular  range  up  to  the  truncation  angles.  This  was 
the  procedure  used  in  the  "end-point  corrected"  seismogram  in  Figure  5.4.  In  the  slab 
problem,  truncation  phases  can  be  easily  reduced  by  including  more  rays  in  a  region  of 
the  model  in  which  the  calculation  can  be  peformed  analytically  and  cheaply.  Cormier 
(1986)  describes  how  fast,  analytic  integration  of  the  dynamic  ray  tracing  equations  can 
be  performed  in  1-D,  radially  symmetric  portions  of  the  modei  and  the  and  simply 
connected  with  numerically  integrated  quantities  in  the  laterally  varying  portions  of  the 
model,  radially  symmetric  portions  in  the  model  is  described  in  Cormier  (1986a). 

The  example  calculations  shown  in  this  proposal  have  not  yet  included  source 
radiation,  attenuation,  reflection/transmission  coefTicents,  and  S  wave  polarization.  All 
of  these  calculations  are  well  developed  and  described  in  a  recent  text  by  Cerven^ 
(1985b).  The  programs  with  which  example  calculations  shown  in  this  proposal  have 
been  made  incorporate  all  of  these  effects,  but  they  have  not  yet  been  all  tested. 
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